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ABSTRACT 


With  recent  delivery  of  littoral  combat  ships,  the  impact  of  operating  in  shallow  or  littoral 
ocean  domain  (LOD)  during  the  duration  of  their  life  cycle  is  of  interest,  and  a  shock  trial 
or  hardening  test  and  validation  for  this  class  is  needed.  For  this  study,  the  theories  of 
underwater  shock  phenomena  as  applied  within  the  boundaries  of  LOD,  specific  to  the 
Eulerian  fluid  domain  were  conducted.  The  results  of  varying  ocean  depth  show  clear 
distinction  in  UNDEX  characterization  at  depths  shallower  than  300ft.  Varying  charge 
size  and  depth  showed  that  charge  size  of  less  than  3001bs  of  HBX-1  displayed  a  linear 
relationship  while  changing  the  charge  depths  to  near  water-air  or  water-bottom  interface 
also  resulted  in  amplified  characteristics  of  UNDEX  parameters.  In  addition,  varying 
lateral  boundary  showed  that  as  its  distance  is  brought  inside  the  radius  of  bulk 
cavitation,  the  UNDEX  behavior  also  became  increasingly  chaotic  due  to  similar  effects 
seen  in  the  shallower  bottom  depth.  Eastly,  adding  blocked  cells  prior  to  a  full  scale 
coupled  run  showed  that  fluid  behaves  more  erratically  as  these  small  rigid  boundaries 
are  situated  within  the  radius  of  forming  bulk  cavitation. 
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I.  INTRODUCTION 


A,  BACKGROUND 

The  theory  of  underwater  explosion  (UNDEX)  phenomena  has  been  studied 
extensively  sinee  post  World  War  II.  For  United  States  Navy  (USN),  its  motivation  was 
to  analyze  damages  incurred  on  surface  and  subsurface  vessels  during  direct  impact  and 
near  misses  from  underwater  weapons.  Since  then,  advancement  in  computing  power, 
extensive  understanding  of  the  UNDEX  theory  and  theory  of  finite  element  have  led  to 
solidifying  efforts  in  validating  shock  trial  analysis  of  newly  commissioned  ships  via 
modeling  and  simulation  prior  to  an  actual  at  sea  shock  trial  analysis.  The  result  is 
software  that  utilizes  UNDEX  theory,  based  on  theoretical  works  of  Arons,  Cole,  Snay 
and  Taylor,  called  Dynamic  System  Mechanics  Advanced  Simulation  (DYSMAS)  [1]. 
DYSMAS  is  a  result  of  the  combined  efforts  of  German  Defense  Contractor, 
Industrieanlagen-Betriebsgesellschaft  mbH  (lABG)  and  Naval  Surface  Warfare  Center 
Indian  Head  Division  (NSWC  -  IHD),  as  well  as,  Eawrence-Eivermore,  led  by  the  USN. 
It  is  a  six  degree  of  freedom  (DOF)  finite  element  hydrocode  used  to  model  and  simulate 
fully  coupled,  fluid-structure  interaction  (ESA)  problems  of  an  UNDEX  event  on  a 
surface  or  subsurface  vessel. 

With  the  recent  delivery  of  the  littoral  combat  ship  (ECS),  the  impacts  of 
operating  in  a  shallow  or  littoral  ocean  domain  (EOD)  during  the  duration  of  the  ship’s 
life  cycle  are  of  interest.  Thus,  a  shock  trial  or  hardening  test  and  validation  for  this  class 
is  needed.  However,  due  to  current  Department  of  Defense’s  (DOD)  budgetary 
constraints  and  rising  ECS  class  issues,  an  actual  shock  trial  analysis  for  ECS  class  has 
been  postponed  indefinitely  [2].  Nevertheless,  continued  study  of  ESA  within  the  EOD 
during  UNDEX  is  needed  in  order  to  verify  the  sea  and  battle  worthiness  of  ECS  class. 

From  Naval  Postgraduate  School’s  (NPS)  Shock  and  Vibration  Computation  Eab 
(SVCE),  several  studies  have  already  been  conducted  to  address  the  factors  of  UNDEX 
simulations  of  littoral  or  shallow  water  by  Walters,  Arbogast,  and  Santiago  regarding 
explicitly  modeled  solid  Eagrangian  ocean  floor,  influence  of  shock-induced  air  bubble 
collapse  and  cavitation  effect  on  surface  ship  model,  respectively.  As  Walters  suggests. 
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the  distinct  advantages  of  using  Lagrangian  solid  bottom  was  the  ability  to  model 
contoured  bottom  shapes  and  for  the  various  solid  bottom  types  created,  there  were 
significant  differences  between  the  initial  bottom  reflections  for  the  different  contours 
[1].  He  further  suggests  that  the  most  important  bottom  contour  effect  was  the  distortion 
to  the  gas  bubble  and  its  associated  first  pulse  timing.  Santiago’s  study  of  bulk  cavitation 
(BC)  effect  on  fluid-structure  interaction  (FSA)  shows  that  the  effects  of  cavitation 
closure  can  be  quite  significant  and  inflict  another  impulsive  force  upon  the  ship  as 
reported  empirically  during  an  actual  shock  trial  report  of  DDG  53,  in  addition  to  the 
incident  wave,  bubble  pulsations,  bottom  reflection,  or  bottom  refractions  [3].  Lastly, 
Arbogast’s  work  showed  that  the  bubbles  in  the  vicinity  of  marine  vessels,  depending  on 
the  size  and  location  and  proximity  of  these  bubbles,  created  a  buffering  effect  from 
shock  induced  by  an  UNDEX  [4], 

B,  SCOPE  OF  RESEARCH 

The  importance  of  FSA  and  depth  of  such  UNDEX  theory  in  evaluating  and 
analyzing  shock  design  of  a  naval  vessel  cannot  be  understated.  However,  in  order  to 
understand  the  FSA,  first,  a  solid  understanding  of  how  fluid  behaves  during  UNDEX 
must  be  clearly  sought  out.  As  Figure  1  shows,  there  are  numerous  theoretical  and 
experimental  analysis  that  deals  with  the  interdisciplinary  studies  of  Naval  Ship  Shock 
Design  and  Analysis  [5].  In  short,  it  is  primarily  divided  into  structural  analysis  that 
stems  from  a  thorough  understanding  of  dynamic  fluid  characteristics  and  fluid-structure 
interactions  during  UNDEX  through  the  use  of  finite  element  modeling  and  simulations, 
as  well  as,  real  world  experimental  validations.  For  this  study,  the  theories  of  underwater 
shock  phenomena  as  applied  within  the  boundary  of  shallow  water  or  EOD  will  be  of 
focus. 

Building  on  the  knowledge  of  previous  researches  conducted  here  at  NFS  and 
various  UNDEX  studies  published  over  the  years,  a  more  in-depth  analysis  and 
characterization  of  EOD  is  attempted.  Furthermore,  a  clear  definition  of  EOD,  in  terms 
of  UNDEX  parameters,  especially  the  influence  on  BC  is  sought  out.  In  addition,  the 
implementation  of  various  boundary  conditions  for  both  bottom  and  side  conditions  of 
fluid  domains  were  compared  and  contrasted  in  the  current  Eulerian  fluid  model,  gaining 
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better  understanding  of  their  effect  on  shallow  water  UNDEX  phenomena.  In  short, 
several  variations  of  typical  benchmark  LOD  problems  were  modeled  and  simulated 
using  DYSMAS  to  gain  better  insights  to  its  resulting  UNDEX  parameters  including  the 
impacts  that  BC  have  on  the  overall  response  to  fluid  behaviors. 
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II.  UNDERWATER  EXPLOSION  THEORY 


A.  INTRODUCTION  TO  CONVENTIONAL,  NON-NUCLEAR  EXPLOSION 

A  conventional,  non-nuclear  explosion  is  a  chemical  reaction  in  a  material  or 
substanee  that  converts  the  original  material  into  high  pressure  gas  and  temperature  at 
extreme  and  catastrophic  speed  [6].  In  general,  any  explosive  material,  regardless  of  its 
state,  is  an  unstable  compound  that  becomes  more  stable  during  the  chemical  reaction 
that  takes  place  post  ignition  and  detonation.  The  temperature  within  the  explosion  of  a 
gas  is  usually  in  the  order  of  3,000°  C  (5,432°  F)  and  the  pressures  are  greater  than 
50,000atm  (734,797psi)  [6].  A  reaction  of  this  magnitude  can  only  be  initiated  with 
sufficient  energy  provided  to  a  specific  place  within  the  explosive  prior  to  detonation  and 
is  done  usually  by  means  of  a  fuse. 


Figure  2.  Conventional,  Non-nuclear  Explosion  on  Land,  from  [7] 

During  detonation  process,  the  developing  intense  heat  and  pressure  are  sufficient 
to  set  up  the  explosive  material  surrounding  the  fuse.  As  a  result,  a  very  narrow 
boundary  between  the  exploding  material  and  its  surrounding  develops  as  a  product  of 
extremely  high  temperature,  pressure  and  energy.  This  distinctly  defined  rapidly 
advancing  discontinuity,  known  as  a  “detonation”  or  “blast  wave,”  travels  with  a  velocity 
of  several  thousand  feet  per  seeond  and  destroys  everything  in  its  path  (Figure  2)  [6]. 
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The  way  in  which  blast  waves  or  detonating  disturbance  propagates  depends  on  the 
physical  or  chemical  properties  of  explosive  materials  and  the  surrounding  medium — air, 
water  or  vacuum — at  which  the  explosion  takes  place.  From  the  standpoint  of 
destructiveness  of  explosive  materials  and  impact  to  its  surrounding,  the  process  of 
detonation,  propagating  shock  waves  (SW)  and  energy  transfer  that  takes  place  through 
the  means  of  pressure  or  velocity  changes  to  its  surrounding  becomes  increasingly 
important. 

B,  UNDERWATER  EXPLOSION 

Conventional  UNDEX  is  a  non-nuclear  explosion  that  takes  place  at  the  surface  or 
subsurface  of  water.  The  characteristics  of  UNDEX  is  similar  to  that  of  conventional 
explosion  discussed  previously,  except  due  to  water  as  medium,  the  initial  mass  of 
explosive  becomes  unstable  hot  mass  of  gas  at  tremendous  pressure  that  radiates 
spherically  affecting  its  surrounding  (see  Eigure  3(a)).  The  energy  expenditure  during 
UNDEX  detonation  accounts  for  predominantly  the  SW  energy  (45%)  and  bubbles  (45%) 
with  approximately  10%  unaccounted  for  (see  Eigure  3(b))  [8].  Within  the  energy 
expenditure  partition  of  SW  and  bubbles,  local  and  BC  are  also  accounted  for  some  of  the 
energy  that  was  imparted  into  its  surrounding  water  during  detonation. 

Eor  this  research,  water  is  characterized  as  a  homogeneous  fluid  capable  of 
supporting  only  limited  shearing  stresses;  hence,  a  medium  and  its  volume  of  which  can 
readjust  itself  to  displacement  of  boundaries  by  result  of  flow  [6].  In  addition,  as 
conservation  of  energy  and  work  implies  the  change  in  its  pressure  on  confined  mass 
results  in  compression  or  work  done  to  that  system,  which  creates  change  in  volume  of 
such  mass.  Hence,  when  pressure  is  applied  to  a  localized  region,  the  compressibility  of 
water  makes  it  possible  for  it  to  transmit  energy  as  a  wave  disturbance  to  other  parts 
within  the  domain  with  a  finite  velocity  which  results  in  a  local  motion  of  water  and 
pressure  changes. 
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Figure  3.  (a)  Shock  Wave  and  (b)  Energy  Partition  at  Time  of  Detonation, 

from  [5],  after  [8] 

If  the  pressure  applied  to  localized  regions  within  the  water  domain  is  small 
enough,  the  rate  of  propagation  or  wave  of  disturbances  can  be  independent  of  its 
pressure  magnitude.  This  propagating  SW  travels  at  about  4,900ft/sec  (fps),  in  sea  water 
at  18°  C  (64.4°  F)  [6].  Of  note,  this  velocity  of  the  propagating  disturbance  is  also 
dependent  on  other  characteristics  of  the  water,  like  salinity  and  its  density.  If  this  wave 
motion  is  simplified  to  one-dimension  so  that  plane  waves  are  assumed,  vice  spherical, 
the  wave  travels  without  significant  changes  in  its  characteristics.  However,  if  the 
detonating  waves  are  radiated  from  a  spherical  source,  there  is  an  inverse  relationship 
between  the  amplitude  of  the  radiating  SW  and  distance  from  its  source.  During  the 
diverging  propagation  of  spherical  wave,  motions  of  its  surrounding  water  are  affected 
due  to  “surge  or  afterflow”  that  takes  place  in  reaction  to  propagating  waves  [6]. 

As  this  SW  propagates,  the  effects  it  has  on  the  surrounding  medium  and  various 
boundaries  like  water-air  and  water-bottom  interfaces,  cause  additional  tensile  or 
compressive  wave  that  induces  additional  havoc  to  existing  underwater  objects.  When 
initial  SW  reaches  water-air  interface,  the  differences  in  tensile  and  compressive 
strengths  of  water  and  air  causes  massive  cavitation  called,  BC  as  it  spallates  the  water 
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particles  on  the  surfaee  and  separates  this  layer  from  below.  The  initial  bubble  with  its 
detonating  eycles  of  its  expansion-collapse  motion  as  it  rises  to  the  surface  and  spews  out 
exploding  material  into  the  atmosphere  due  to  buoyaney,  it  too  creates  another 
catastrophie  havoe  to  anything  in  its  path  or  in  the  vicinity  there  of.  Hence,  the  five 
primary  events  that  take  place  upon  detonation  during  UNDEX  are  initial  SW 
propagation,  formation  of  BC,  cyclic  bubble  phenomena,  seeondary  or  third  pressure 
pulses  due  to  boundary  reflections  and  surface  phenomenon. 

1.  The  Shock  Wave  and  Acoustic  Wave  Assumptions 

The  initial  cause  of  disturbance  to  the  water  in  an  UNDEX  is  the  arrival  of  shook 
pressure  wave  reacting  between  the  boundaries  of  explosives  and  water  [6].  This 
pressure  distribution,  usually  in  the  order  of  34,0001b/in  (psi)  for  a  3001bs  of  TNT, 
begins  to  rapidly  decay  as  it  propagates  followed  by  intense  peak  pressure  wave  in  an 
outward  motion  of  the  water  (figures  4  and  5)  as  the  extremely  dense  mass  of  exploding 
gas  oontinues  to  spherioally  expand  as  detonation  takes  places.  As  its  SW  pressure 
diminishes,  so  too  does  the  pressure  in  the  water  fall  off  rapidly,  as  shown  in  three 
different  target  distances  of  5,  50  and  500ft  from  the  point  of  detonation  (Eigure  4). 
Easting  for  a  few  milliseconds  at  most,  a  discontinuity  exists  in  the  initial  SW  pressure 
rise  that  is  followed  by  a  distinct  exponential  decay,  shown  in  Eigure  5  [6]. 


Eigure  4.  Pressure  Distribution  at  (a)  5  ft  (b)  50ft  and  (c)  500ft  from 

Detonation,  from  [6] 
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The  pressure  distributions  of  a  3001bs  TNT  charge,  shortly  after  detonation  at 
target  distance  of  50  and  500ft  at  three  instants  of  time  from  the  center  of  TNT  charge  are 
depicted  in  Figure  5.  For  comparison,  the  initial  wave  in  previous  Figure  4(a)  is  also 
depicted  in  acoustic  wave  form  by  dotted  lines  for  (b)  50  and  (c)  500ft.  The  obvious 
scaling  law  that  emerges  when  different  charge  sizes  are  used  is  called,  “principle  of 
similarity.”  The  principle  of  similarity  states  that  “if  the  linear  size  of  charge  is  changed 
by  factor  of  k,  the  pressure  conditions  will  also  change  at  new  distance  and  time  scales  k 
times  as  large  as  the  original  ones  are  used”  [6]. 


Figure  5.  Shock  Wave  Pressure  Profile  at  Two  Distances  from  Detonation, 

from  [6] 


The  radially  propagating  disturbance  as  a  wave  of  compression  in  water  called, 
“shock  wave”  has  following  characteristics  [6]; 

a.  The  initial  velocity  of  propagation  at  or  near  charge  and  water  boundary  is 
approximately  5,000fps  at  detonation  and  drastically  falls  to  its  “acoustic” 
values  as  wave  advances  outward  (Figure  4) 

b.  The  pressure  level  in  spherically  propagating  wave  falls  off  more  rapidly 
than  with  inverse  of  distance  near  the  charge  at  detonation,  but  eventually 
converges  to  this  behavior  as  it  reaches  larges  distances. 

c.  The  region  of  highest  peak  pressure  is  at  target  location  closest  to  the 
charge  at  detonation  and  its  wave  profile  broadens  gradually  as  the  wave 
propagates  outward.  Same  characteristics  can  be  observed  at  various  target 
location  with  reduced  peak  pressure. 
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As  mentioned,  sinee  the  initial  SW  veloeity  will  be  drastieally  redueed  and 
eventually  eonverged  to  aeoustie  speed  as  it  propagates,  it  is  assumed  that  the  SWs  are 
approximated  as  aeoustie  plane  waves  with  small  compression  and  speed.  As  a  result, 
following  Equation  (2.1)  is  derived. 


(2.1) 


which  states  that  the  added  fluid  pressure  (P)  from  propagating  shockwaves  is  a  function 
of  water’s  mass  density  (po),  acoustic  wave  speed  (Co)  and  its  propagating  velocity  (u). 
The  medium  in  which  SW  travels,  water  in  this  case,  are  considered  to  be  constant 
properties  and  as  mentioned  previously,  this  relationship  assumes  that  water  is  uniform, 
homogenous,  compressible  and  inviscid. 


Figure  6.  UNDEX  Geometry  and  Various  Paths  of  Detonating  Shockwave, 

after  [5] 

There  are  basically  four  primary  paths  that  a  propagating  shockwave  can  travel 
during  an  UNDEX  event  before  reaching  its  target  (Figure  6):  (1)  Direct  Path  (2)  Surface 
Reflected  or  Rarefacted  (3)  Bottom  Bounce  and  (4)  Bottom  Refracted  Waves. 
Depending  on  the  UNDEX  geometry  of  the  explosion  in  terms  of  charge  depth,  size  and 
ocean  depth,  these  shockwave  travels  in  the  order  of  listed  and  the  nature  of  applied 
pressure  waves  are  either  compressive  (+)  or  tension  (-)  as  experienced  by  the  target  due 
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to  varying  densities  at  water,  air  or  bottom  interfaees.  Of  note,  the  tensile  pressure  wave 
of  (2)  surfaee  eut-off  is  due  to  water’s  inability  withstand  tensile  foree  at  air-water 
interfaee;  henee,  ereating  a  region  of  expansion  in  water  moleeule  rather  than  a 
eompressive  stress  near  detonation  and  other  traveling  path  of  pressure  wave. 


Figure  7.  Arbitrary  Interfaee  Geometry,  from  [1] 

The  initial  shoekwave’s  behavior  at  an  interfaee  sueh  as  water- air  or  water-bottom 
(seabed)  ean  be  approximated  by  using  eombinations  of  Snell’s  law  and  pressure-veloeity 
eontinuity  at  the  interfaee  [5].  An  arbitrary  two  different  medium  interfaee  is  shown  in 
Figure  7.  Using  Snell’s  law,  the  relationship  of  various  medium’s  aeoustie  speed  (C)  and 
its  angle  of  ineident  (a)  to  the  approaehing  interfaee  ean  be  refereneed  to  its  normal 
veetor  of  that  interfaee  or  summarized  in  Equation  2.2  below: 


C, 


sin  n 


sin 


(2.2) 


Now,  if  the  pressure-veloeity  eontinuity  is  applied,  equations  (2.1)  and  (2.2)  ean 
be  eombined  to  deseribe  the  interfaee,  yielding  the  following  relationship  between  the 
pressure  and  ineidental  veloeity,  as  well  as,  refleeted,  refraeted  or  transmitted 
shoekwaves. 


11 


p 

^refl 

P. 

me 


P2C2  cos  cos  ^trans 

P2C2  COS  cos  ^trans 


U. 


‘refl 


Uinc 


(2.3) 


H, 


2^2(22  cos 


P, 

prjc 


^2(22  COS  cos 


P2^2  Pt 
u 


(2.4) 


Therefore,  when  above  relationships  are  applied  to  varying  boundaries,  the 
resulting  pressure  wave  beeomes  either  eompressive  or  tensile  waves.  At  rigid  boundary 
(P2C2  »  PiCi),  sueh  as  the  water-bottom  interfaee,  the  wave  refleetions  are  compressive 
with  the  magnitude  of  both  incident  and  reflected  waves  pressure  and  velocity  being 
equal,  since  by  definition,  the  velocity  of  a  rigid  boundary  is  zero  [1].  For  a  non-rigid 
boundary  condition  (P2C2  «  piCi),  such  as  water-air  interface,  the  same  relationships 
applied  as  above  results  in  a  tensile  reflection  or  rarefaction  waves. 
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Figure  8.  MATLAB  Figure  of  Pressure  History,  from  [1] 


The  simultaneous  passage  of  a  compressive  and  tensile  wave  through  a  point  at 
different  times  causes  a  unique  effect  creating  a  discontinuous  drop,  is  called  the  pressure 
cutoff  shown  in  pressure-time  history  in  Figure  8.  At  this  point,  the  abrupt  halt  to  the 
impulse  induced  on  a  target  at  the  point  of  cut-off  can  also  be  observed  [1].  This 
dependence  of  UNDEX  behaviors  on  its  geometry  is  depicted  in  Figure  9  and  using  such 
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geometry,  pressure  and  vertical  take-off  velocity  (VTO)  can  be  calculated  to  analyze 
target  behavior  at  various  instants  of  time. 


Figure  9.  UNDEX  Geometry,  from  [1] 

The  obvious  interdependence  of  charge  weight  and  wave  propagation  distance  for 
spherical  charge  gave  rise  to  empirical  formulas  to  approximate  several  important 
parameters  of  UNDEX  that  have  been  continuously  developed  and  tested  by  Cole  and 
others  since  WWII  [5].  The  two  primary  formulas  are  the  peak  pressure  (Pmax)  and  decay 
constants  (0)  shown  in  Equations  (2.5)  and  (2.6): 

(2.5) 

(2.6) 

where  charge  weight  (W)  is  in  pounds,  the  wave  propagation  distance  (R)  is  in  feet,  the 
peak  pressure  (Pmax)  in  psi  and  the  decay  constant  (0)  in  milliseconds  (msec).  The 
constants  K  and  A  are  depended  on  the  explosive  materials  used  and  were  found 
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experimentally  prior  to  this  researeh  through  numerous  field  testing  during  various 
historieal  researehes  [5]. 

2.  Typical  UNDEX  Problem 

A  typieal  UNDEX  problem  with  an  explosive  eharge  of  weight,  W,  detonated  in 
water  at  depth,  D,  is  eonsidered  in  Figure  10  and  its  eorresponding  pressure  wave  profile 
as  time  passes  is  shown  in  Figure  11  [9].  The  water  depth  for  this  problem  is  assumed  to 
be  infinitely  deep;  therefore,  negligible  when  it  eomes  to  bottom  bounee  interaetion  and 
eonsiderations.  In  order  to  measure  various  parameters  of  UNDEX  event  taking  plaee,  an 
absolute  pressure  gage  is  loeated  at  target  loeation  of  point  (X,  Y)  or  1  to  measure  the 
pressure  ehanges  at  this  loeation.  As  detonation  takes  plaee,  a  eompressive  SW 
propagates  radially  or  spherieally  from  its  souree  and  travels  along  the  line  0  -  1  in 
Figure  10  as  it  reaehes  the  target  or  gage  loeation  [9]. 


Figure  10.  Ineident  and  Refleeted  Wave  Paths,  from  [9] 

The  absolute  pressure  gage  shows  eonstant  pressure  at  its  hydrostatie  depth 
denoted  by  line  A  -  B  in  Figure  1 1  prior  to  detonation.  This  pressure  at  this  depth  is  sum 
of  atmospherie  and  hydrostatie  pressure  at  the  gage  depth  and  upon  arrival  of  SW  at 
target  point,  the  pressure  jumps  to  point  C,  followed  by  exponential  deeay  along  line  C  - 
D  is  observed  [9].  As  will  be  shown  for  most  of  the  eases  studied  in  this  researeh,  this 
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pressure  is  at  or  close  to  14.7  psi  or  latm;  for  all  practical  purpose,  at  or  near  surface  of 
the  water.  As  shown  previously  in  Figure  6  and  now  Figure  10,  the  line  0-2  represents 
the  compressive  path  of  SW  as  it  travels  to  the  water-air  interface.  Next,  the  SW  takes 
the  path  of  2  -  1,  where  a  rarefaction  or  tensile  wave  is  observed  due  to  difference  in 
acoustic  impedance  of  air  and  water  mentioned  previously. 


Figure  1 1 .  Absolute  Pressure  at  a  Point,  from  [9] 


The  arrival  of  this  rarefaction  or  surface  reflected  wave  at  the  gage  location  is 
termed  surface  cutoff  and  is  illustrated  in  Figure  1 1  by  a  sudden  drop  in  the  absolute 
pressure  from  point  D  to  three  different  locations,  E,  F  or  G,  depending  on  the  vapor 
pressure  of  the  surrounding  water  [9].  This  concept  of  surface  cutoff  is  made  clearer  by 
depicting  the  rarefaction  wave  as  emanating  from  an  image  charge,  Wi,  and  propagating 
along  line  3  -  2  -  1  in  Figure  10.  This  method  of  including  image  charge  to  calculate 
surface  cutoff  characteristics  is  known  as  method  of  images  and  is  valid  since  the 
distance  3-2-1  equals  the  distance  0  -  2  -  1 .  Therefore,  where  point  D  drops  upon 
cutoff  as  the  rarefaction  wave  tries  to  lower  the  absolute  pressure  at  the  gage  location  is 
depended  on  one  of  three  levels  [9]; 

E:  If  pressure  at  E  is  greater  than  vapor  pressure 

E:  If  pressure  at  E  is  great  than  cavitation  but  less  than  vapor  pressure 

G:  If  pressure  at  G  is  less  than  cavitation  pressure 
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3. 


Bulk  Cavitation 


The  arrival  of  rarefacted  or  reflected  wave  near  surface,  also  called  cut-off  due  to 
water’s  inability  to  withstand  tensile  force  at  air-water  interface,  creates  BC  near  air- 
water  interface.  In  general,  unlike  seawater,  pure  homogenous  water  with  no  impurities 
is  able  to  withstand  a  considerable  amount  of  tensile  stress.  Due  to  abundance  of  obvious 
biological  and  non-biological  impurities,  seawater  can  withstand  very  little  tension.  As  a 
result,  the  layer  of  sea  water  at  or  near  the  water-air  interface  ruptures  and  the  pressure 
falls  to  the  vapor  pressure  of  water,  assumed  to  be  approximately  zero  for  all  practical 
purpose  when  dealing  with  UNDEX  [10].  Within  the  areas  of  rupturing  water  layer  near 
the  surface,  a  collective  spallation  takes  places  vertically  and  tensile  stresses  at  its  lower 
boundary  tending  the  water  particles  to  return  to  its  original,  pre-shocked  position  cease 
to  exist.  Of  note,  the  motion  of  spallation  is  dictated  primarily  by  the  atmospheric, 
hydrostatic  pressure  and  the  acceleration  due  to  the  gravity.  The  radial  and  vertical 
thickness  of  typical  BC  zone  for  a  601bs  charge  of  HBX-1  at  a  depth  of  30ft  is  depicted  in 
Figure  12.  As  shown,  this  charge  size  and  depth  creates  a  BC  that  is  574ft  in  radius  and 
24.6ft  deep  at  its  maximum  thickness;  hence,  as  shown,  its  radius  is  several  magnitudes 
longer  than  its  radius. 


Figure  12.  BC  Zone  for  60  lbs  of  HBX-1  at  30ft 


Next,  the  spalled  water  layer  is  driven  back  downward  by  force  of  gravity  so  that 
ultimately  it  must  impact  the  underlying  water  and  produces  a  traveling  source  of 
secondary  pressure  wave  in  the  water  called  “hammer.”  This  traveling  source  and 
hammer  created  are  almost  equivalent  to  a  “sonic  boom”  in  the  atmosphere  and  because 
of  its  focusing  effects,  the  secondary  pressure  due  to  hammer  can  reach  intensity  as  high 
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as  initially  propagating  SW  [10].  Interestingly,  the  region  of  cavitation  between  the 
upward  traveling  spalled  water  layer  and  the  relatively  quiescent  underlying  body  of 
water  have  an  average  specific  volume  greater  than  that  of  seawater  and  this  intervening 
region  is  termed  the  BC  region  due  to  its  size  and  can  be  observed  in  figures  12  and  13. 


Figure  13.  Photograph  of  Upper  Surface  of  Bulk  Cavitation,  from  [11] 


Finally,  in  an  attempt  to  summarize  numerous  BC  studies  conducted  and  theories 
derived  over  the  years  since  the  1960s,  the  critical  sections  within  the  technical  paper 
published  by  Costanzo  and  Gordon  in  May  1983,  titled  “A  Solution  to  the  Axisymmetric 
Bulk  Cavitation  Problem”  [9],  as  well  as,  sample  case  studies  are  all  included  in  the 
Appendix  of  this  study  as  a  reference  and  also  to  be  used  in  subsequent  chapter’s  analysis 
of  LOD  UNDEX  phenomenon  for  its  underlying  characteristics  of  BCs. 

4,  Motion  of  Gas  Sphere  and  Its  Empirical  Bubble  Formula 

During  the  expansion  and  collapse  of  cyclic  bubble,  the  initial  peak  pressure  in 
the  expanding  gas  sphere  is  decreased  dramatically  due  to  energy  dissipation  to  its 
surrounding  from  detonation.  However,  the  immediate  state  of  the  spherical  bubble  that 
follows  still  has  much  higher  pressure  than  the  equilibrium  hydrostatic  pressure  at  that 
depth  [6].  This  sphere  or  bubble  has  a  tremendous  outward  velocity  as  its  diameter 
increases  rapidly  and  is  proportional  to  the  pressure  magnitude  at  that  exact  time.  The 
expansion  of  bubble  continues  for  a  relatively  long  period  of  time  compared  to  other 
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transient  UNDEX  phenomenon  like  initial  shoek  waves  or  BC.  As  this  bubble  reaehes  its 
maximum  radius,  the  internal  gas  pressure  deereases  dramatieally,  while  the  expanding 
motion  eontinues  due  to  the  inertia  of  water  flowing  outward  [6].  As  the  internal  gas 
pressure  eventually  falls  below  the  equilibrium  pressure  (hydrostatie  +  atmospherie),  this 
outward  motion  or  flow  stops  and  the  bubble  begins  to  contraet  at  an  increasing  rate.  The 
contraction  or  inward  motion  of  the  bubble  continues  until  it  reaches  the  maximum 
compressibility  of  the  gas  and  reverses  the  motion  once  again  due  to  reloading  of 
pressure  differences  within  the  bubble  and  its  surrounding.  Hence,  the  created  elastic 
properties  of  the  gas  and  water  due  to  detonating  pressure  and  hydrostatic  influences  of 
surrounding  water,  in  time  creates  a  condition  of  oscillating  system  that  repeats  itself 
until  it  reaches  the  surface  of  the  water  [6].  This  cyclic  nature  of  bubble  expansion  and 
collapse  as  function  of  time,  as  well  as,  its  displacement  as  it  travels  up  towards  the 
surface  due  to  hydrostatic  forces  induced  by  surrounding  water  and  gravity  is  depicted  in 
Figure  14. 


Figure  14.  (a)  Pulsation  and  (b)  Displacement  of  Gas  Bubble,  from  [6] 


The  traveling  cyclic  bubble  pulses  and  its  comparison  to  corresponding  pressure 
profile  during  FfNDEX  are  shown  in  Figure  15.  As  mentioned,  the  initial  shockwave 
pressure  is  the  greatest  at  minimum  bubble  radius  at  the  time  of  detonation.  As  the  cyclic 
bubble  pulse  travels  to  near  water-air  interface,  the  detonating  energies  are  dissipated  and 
the  bubble  eventually  reaches  the  surface  to  create  plume  into  the  atmosphere.  This 
nature  of  cyclic  bubble  can  induce  hull  whipping  to  surface  vessels  if  period  of  the 
bubble  pulse  matches  its  lowest  natural  frequency;  hence,  the  importance  of  analyzing 
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such  phenomenon  during  UNDEX.  Within  the  energy  partition  of  the  shockwave  bubble, 
the  initial  shoekwave  dominates  initial  energy  dissipation  followed  by  cyelie  bubble 
pulses  (Figure  16).  The  initial  SW  energy  can  be  as  high  as  50%  of  total  bubble 
dissipated  energy  while  the  remaining  energies  are  attributed  to  traveling  eyclic  bubbles 
pulses  [8]. 


Figure  15.  Traveling  Cyelie  Bubbles  and  Its  Pressure  Profile,  from  [5] 


Energy  Partition  of  Shockwave  Bubble 

■  Initial  Shockwave 

■  1st  Pulse 

■  2nd  Pulse 


Figure  16.  Energy  Partition  of  Shockwave  Bubble,  after  [8] 

Numerous  theoretical  methods  and  approaehes  to  modeling  and  analysis  of  gas 
bubbles  and  its  oseillations  have  been  published  over  the  years.  However,  for  this  study, 
a  shallow  water  assumption  was  made  to  limit  the  seope  of  sueh  analysis.  As  a  result, 
two  parameters  of  gas  bubble  emerge  during  FTNDEX,  its  maximum  radius  and  period  of 
initial  pulse  [1].  Taking  into  aecount  the  correetion  faetor  that  must  be  ineluded  due  to 
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proximity  of  the  air-water  or  water-bottom  interfaces  in  shallow  water,  as  well  as, 
following  the  same  intuitions  in  deriving  shockwave  pressure  equations  (2.5)  and  (2.6), 
the  maximum  bubble  radius  (Amax)  in  feet  is  derived  [5]: 


^max  -^6 


I¥ 


D  +  33 


(2.7) 


Next,  the  period  of  initial  bubble  pulse  (T)  is  derived  in  seconds  with  its  correction  factor 
due  to  distance  from  surface,  shown  in  Equation  (2.8).  The  correction  factor  of 
numerical  constant  (a)  is  equal  to  0.1  when  Amax/D  is  less  than  0.5  and  when  greater 
than  0.5,  then  a  is  approximated  between  0.1  and  0.2  [8]. 
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5,  Secondary  or  Additional  Pressure  Pulses  and  Surface  Effects 

As  described  previously,  the  bubble  pulse  due  to  cyclic  expansion  of  gas  as  it 
travels  to  the  surface  of  the  water  depends  considerable  on  the  water  depth  and  proximity 
of  the  boundary  surfaces  at  detonation.  The  first  peak  pressure  of  bubble  pulse  is  usually 
no  more  than  10  to  20%  of  that  of  initial  SW  [6].  However,  its  duration  is  much  longer 
and  the  areas  under  the  pressure  profile  curve  are  considerably  close  [6].  During 
completion  of  each  cyclic  bubble  pulse,  a  considerable  amount  of  the  energy  initially 
present  at  detonation  is  lost  due  to  energy  transfer  that  occurs  to  its  surrounding  water  as 
it  travels  and  pulsate,  as  seen  previously  in  figures  14  and  15.  As  a  result,  only  the  first 
bubble  pulse  is  of  practical  significance  since  the  successive  pulses  progressively  weaken 
to  negligent  magnitude.  The  relations  between  initial  SW  and  the  first  bubble  pulse 
pressures  and  durations  are  shown  in  Figure  17  and  Figure  18  shows  close  up  details  of 
the  bubble  pulse  pressure  generated  from  300  pounds  of  TNT  charges  detonated  at 
various  depths  in  100  feet  total  depth  of  water. 
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Figure  17.  Pressure  60ft  from  3001bs  TNT  Charge  Detonated  at  50ft,  from  [6] 


Figure  18.  Various  Bubble  Pulses  of  60ft  at  Various  Detonated  Depth,  from  [6] 

It  is  observed  that  the  profile  of  the  first  bubble  pulse  eurve  in  Figure  18  beeomes 
inereasingly  irregular  for  initial  charge  position  close  to  either  surface  or  bottom  [6]. 
Also,  for  shallow  water  or  LOD  conditions,  the  traveling  cyclic  bubbles  and  the  reflecting 
pressure  waves  from  the  surface  or  bottom  give  rise  to  interference;  therefore,  as  would 
be  observed  in  this  study,  the  later  portions  of  observed  pressure  profdes  are  considerably 
different  than  what  would  be  observed  in  an  infinitely  deep  water  domain  [6].  Lastly,  the 
interaction  on  the  water-air  surface  boundaries  as  initial  SW  and  bubble  reaches  the 
surface  is  also  important.  While  the  initial  SW  creates  BC  at  this  interface,  the  spray 
dome  and  subsequent  bubble  pulses  can  generate  multiple  overlapping  plumes  that 
literally  shoot  or  hurl  the  mixture  of  explosive  materials  and  water  vapors  into  the 
atmosphere  creating  additional  havoc  during  UNDEX.  This  surface  phenomenon  which 
marks  the  end  of  major  UNDEX  events  is  shown  in  Eigure  19. 
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(a)  Spray  Dome  (b)  First  Plume  (cl)  Second  Plume 


Coused  by 
Shock  Wave 


Caused  by 
Rrst  Bubble 
Pulse 


Caused  by 
Second  Bubble 
Pulse 


Figure  19.  Surface  Phenomena  (a)  Spray  Dome  or  Spallation  (Bulk  Cavitation) 
Caused  by  Initial  Shock  Wave  (b)  First  Plume  Caused  by  the  First 
Bubble  Pulse  (cl,  c2)  Subsequent  Plume  caused  by  Burping  Bubble 

Pulses,  from  [5] 
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III.  MODELING  AND  SIMULATION  USING  DYSMAS 


A.  DYSMAS 

The  software  paekage  used  for  this  research  is  the  DYSMAS  finite  element 
hydrocode.  DYSMAS  is  a  fully  coupled;  six  degree  of  freedom  (DOF)  finite  element 
(FE)  based  hydrocode  that  is  designed  to  simulate  three-dimensional  (3D),  UNDEX 
events  to  analyze  fluid-structure  interactions  (ESA).  The  foundation  of  DYSMAS’ 
theoretical  concepts  is  based  on  EE  methods  and  like  any  common  EE  software  suite 
currently  available,  it  consists  of  two  major  fluid  and  structural  solvers,  Gemini  and 
Dyna_N(3D),  that  is  interfaced  by  the  Standard  Coupler  Interface  (SCI)  along  with  a  pre- 
and  post-  processor  called,  DYSMAS/P  [12].  The  basic  DYSMAS  suite  architecture  is 
depicted  in  Eigure  20. 


DYSMAS  Coupled  Code  Architecture 


Dyna_N(3D) 

Massively  parallel  explicit 
finite  elentent  solver 


RIgidBody 

Explicit  rigid  bodydynannics 


Gemini 

Massively  parallel  explicit 
compressible  flow  solver 

Replay 

Rapidly  replay  structural  loading 

Float 

Initialize  draft,  trim,  &  roll 

PreStress 

Statically  load  structure 


Eigure  20.  Basic  DYSMAS  Coupled  Code  Architecture,  from  [12] 


The  Gemini  is  an  Eulerian  fluid  solver  and  the  Dyna_N(3D)  is  a  Eagrangian 
structural  solver  that  evaluates  the  structural  dynamic  response  during  UNDEX.  The  SCI 
is  the  key  coupler  that  connects  and  interfaces  during  a  fully  “coupled”  simulation  of  both 
solvers  by  passing  information  between  the  two  solvers  at  the  end  of  each  Eulerian  time 
step,  maintaining  fully  coupled  interface.  There  are  several  other  structural  solvers  that 
can  be  used  in  coupled  DYSMAS  run,  such  as  RigidBody  solver.  DYSMAS/P  and  other 
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post  processors  like  MATLAB  were  used  for  detailed  analysis  on  the  results,  both 
graphieally  and  numerieally. 

1.  GEMINI 

For  this  study,  only  the  Gemini  portion  of  DYSMAS  was  utilized  in  Eulerian  fluid 
behavior  analysis  during  UNDEX.  The  Gemini  eode  is  a  finite-differenee,  Eulerian  fluid 
equation  solver  that  solves  fluid  mesh  by  running  a  higher  order  Godunov  method 
algorithm  to  solve  the  fluid  domain  at  eaeh  time  step  [13].  Designed  speeifically  to 
simulate  the  shockwaves,  cavitation  and  bubble  jetting  phenomena  of  ETNDEX,  Gemini  is 
capable  of  solving  flow  fields  involving  several  Eulerian  fluid  types.  Eising  conservation 
of  mass,  momentum,  and  energy,  the  Godunov  method  algorithm  is  applied  to  solve  eaeh 
Eulerian  equation  for  each  cell  within  the  fluid  domain  with  a  major  assumption  of  Euler 
equation  being  that  the  flow  is  frictionless  or  inviseid. 


Eigure  21 .  Basie  Components  of  GEMINI,  the  Eulerian  Solver,  from  [12] 

Gemini  eonsists  of  several  additional  components  eodes:  GemGrid,  PreGemini, 
Float,  and  Prestress  are  preproeessors.  The  user  defined  and  required  “input”  files 
(shown  in  blue.  Figure  21)  are  used  to  speeify  each  component  eodes.  The  sample  input 
files  are  also  included  in  the  Appendix  of  this  report.  These  sub  component  codes  are 
used  to  eustomize  fluids  domain  set  up  for  a  given  ETNDEX  problem.  GemGrid  is  used 
to  set  up  the  initial  Euler  cell  grid,  up  to  three-dimensions  and  can  be  used  for  any  grid 
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refinement  that  is  needed  around  speeifie  areas  of  interest  in  the  flow  field,  speeifieally 
the  fluid  bubble,  eavitation  and  the  Lagrangian  struetures.  PreGemini  is  used  to  fill  the 
Euler  cell  grid  with  the  user-defined  Eulerian  materials  that  utilizes  the  given  “material 
files”  that  encompasses  specific  material’s  equation  of  state  and  initial  properties,  such  as 
energy,  density,  and  etc.  In  order  to  accurately  position  and  prestress  the  Dyna  N 
structure  in  the  fluid  domain  prior  to  transient  analysis,  the  Eloat  and  Prestress  processors 
are  used. 

As  Gemini  solves  each  Eulerian  cell  grid,  it  dumps  and  saves  all  flow  field  and 
historical  data  into  a  bin  file,  which  in  turn  are  finally  processed  by  GemHis  and 
Gemfield  processor.  GemHis  processes  the  user  specified  recorded  ETNDEX  variables 
and  locations  within  the  fluid  domain  for  analysis.  Gemfield  generates  contour  plots  for 
the  similar  user  defined  ETNDEX  variables  of  designated  plane  snap  shot  at  specified  time 
increments.  Once  flow  field  and  historical  data  files  are  generated  using  GemHis  and 
Gemfield,  DYSMAS/P,  an  independent  DYSMAS  pre/post  processing  suite  or  MATEAB 
can  be  utilized  to  further  analyze  the  Eulerian  and  Eagrangian  solutions.  Samples  of  all 
input  files  are  included  in  the  Appendix  of  this  report. 


B,  DYSMAS/P  AND  MATLAB 

DYSMAS/P  and  MATLAB  were  used  for  pre  and  post  processing  of  DYSMAS 
simulations.  The  preprocessing  of  all  Lagrangian  EE  models  can  be  conducted  using  the 
DYSMAS  Pre-Processor  2010  and  some  Eulerian  characteristics  like  BCs  were  done 
using  MATLAB.  The  preprocessing  capabilities  of  DYSMAS/P  can  be  used  for  the 
creation  of  the  PE  model’s  structure,  assignment  of  material  properties,  boundary 
conditions,  body  forces,  motions,  and  fluid-structure  interface  segment  definitions.  The 
structural  model  is  then  written  into  the  specific  input  cards  required  Dyna  N,  in  turn 
used  as  additional  input  files  during  coupled  runs.  The  DYSMAS/P  postprocessor  allows 
the  user  to  analyze  the  entire  fluid-structure  dynamic  response,  allowing  user  to  extract 
specified  data  and  a  snap  shot  of  specified  time  steps  for  graphical  and  numerical 
analysis.  MATLAB  was  also  used  in  numerical  post  processing  of  DYSMAS  solutions. 
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C.  BASIC  DYSMAS  SIMULATIONS  SETUP 

The  basic  DYSMAS  simulation  set  up  entailed  nearly  the  same  basic  steps  and 
approaches  on  all  simulated  cases  with  minor  differences  depending  on  the  type  of  run  as 
fluid  only  or  fluid-structure  coupled  run.  For  this  study,  fluid  only  runs  were  conducted. 
Following  steps,  in  the  order  as  listed  and  corresponding  input  fdes  were  utilized: 

1 .  Utilize  DYSMAS/P  and  MATLAB  for  preprocessing  as  needed  for 
Lagrangian  FE  structure  model  and  specific  Eulerian  fluids  characteristics 

2.  Generate  Euler  grid  using  Gemgrid 

3.  Eill  the  fluid  domain  using  Pregemini  &  Material  Piles 

4.  Check  fluid  field  prior  to  Gemini  execution  using  Gemfield 

5.  Run  Gemini  (Por  coupled  run,  Gemini  will  be  specified  as  such,  invoking 
fully  coupled  3D  run  using  Dyna_N(3D)) 

6.  Process  Gemini/Dyna_N(3D)  solutions  using  Gemhis  and  Gemfield 

7.  Detailed  post  processing  and  qualitatively  validate  results  prior  to  using 
DYSMAS/P  and  MATLAB  for  final  post  processing 

Like  any  other  PE  methods  or  software,  the  “art”  of  utilizing  PE  was  also 
considered  throughout  the  process,  especially  during  the  mesh  or  grid  set  up,  validation 
of  initial  results  and  distinguishing  the  best  use  of  computational  power  versus  time  for 
each  simulations. 

D,  EULER  GRID,  FLUID  DOMAIN  SETUP  AND  CHARGE  PARAMETERS 

During  grid  setup  using  Gemgrid,  the  most  important  decision  made  were  where 
to  refine  grid  and  where  to  maintain  as  default  coarse  grid  of  approximately  0.656  ft  per 
grid  cell.  According  to  Didoszak  and  from  previous  studies  conducted  at  NPS, 
DYSMAS  simulations  are  highly  sensitive  to  the  mesh  size  surrounding  the  explosion 
[14].  It  was  found  that  the  smaller  the  mesh  or  grid  size,  the  closer  in  accuracy  in  the 
simulation  for  both  empirical  and  experimental  predictions  for  various  UNDEX 
parameters.  Obviously,  these  assumptions  had  to  be  taken  into  consideration  with  a  grain 
of  salt,  because  mesh  refinement  is  not  the  solution  to  all  PE  problems.  In  addition. 
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Gemini’s  free  boundary  conditions  did  not  adequately  allow  fluid  to  flow  in  and  out  of 
fluid  domain  quickly  enough,  if  the  fluid  domain  was  too  small  in  comparison  to  the 
maximum  bubble  diameter.  Lastly,  when  setting  bottom  boundary  conditions,  anything 
less  than  “wall”  or  rigid  gave  unexpected  errors  in  flow  field  and  historical  results.  This 
finding  and  realization  was  particularly  useful  and  had  to  be  considered  in  shallow  water 
simulations  and  solutions  where  overall  fluid  depth  was  restricted  considerably  compared 
to  deep  open  ocean  scenarios. 


Air  Column 

I- - >  X 

\  Charge  (X,  Y,  Z) 

Water  Column 


_ Clay  Column 

Figure  22.  Typical  Fluid  Domain  Setup  for  Gemini 

Following  the  grid  set  up,  the  fluid  domain  was  initialized  and  populated  through 

the  Pregemini  input  file.  As  typical  Eulerian  fluid  domain  setup  is  depicted  in  Figure  22 

and  the  critical  step  in  this  procedure  was  the  seeding  of  the  fluid  domain  with  the  correct 

materials  at  the  appropriate  states  and  depth  or  geometry  and  verifying  it  with  Gemfield 

prior  to  Gemini  execution.  In  order  to  populate  the  fluid  domain  with  different  materials 

such  as  explosives,  air,  water,  etc.  Pregemini  was  utilized  with  given  various  material 

files.  Within  the  material  files  contained  specified  equation  of  state  and  variables 

(pressure,  energy,  and  density)  that  can  also  be  manipulated,  if  needed,  for  specified 

materials  at  user’s  discretion.  As  part  of  the  fluid  domain  set  up  within  the  Pregemini 

input  file,  except  where  specifically  noted,  a  standard  charge  of  HBX-1  was  used.  The 
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charge  weight  and  depth  were  varied  or  kept  eonstant  aeeording  to  the  speeifieity  of  eaeh 
simulation  eonducted. 

E.  SUMMARY  OF  SIMULATION  SETUP 

Taking  all  previously  diseussed  factors  of  UNDEX  theory  and  lesson  learned 
from  initial  simulations  and  previous  studies,  the  deeiding  values,  properties  and  setup  for 
eaeh  simulation  started  from  the  basie  guidance  and  logies  established  in  this  seetion. 
This  basic  set  up  was  then  further  explored  and  developed  into  eomplex  and  more 
realistie  Eulerian  fluid  model  as  eaeh  DYSMAS  solutions  were  qualitative  validated  in 
aceordanee  with  UNDEX  theories  diseussed  in  previous  seetions  of  this  report  prior  to 
final  simulation  runs.  Eaeh  ease  studies  and  simulation  set  up,  including  its  parameters 
and  eontributing  faetors  ean  be  found  throughout  subsequent  ehapters  of  this  report  and 
the  eomplete  outline  of  ease  studies  ean  be  found  in  the  Appendix  of  this  report. 
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IV.  INTRODUCTION  LITTORAL  OCEAN  DOMAIN 


Recent  proliferation  of  asymmetric  threats,  such  as,  diesel  subs,  mines,  fast  patrol 
boats,  coastal  integrated  defense  in  littoral  areas  and  even  drones  makes  littoral  ocean 
domain  (LOD)  a  battle  space  that  USN  must  carefully  consider  and  analyze  in  the  future 
maritime  battle  scenarios  [15].  The  term,  “littoral”  water  (as  depicted  in  Figure  23)  has 
been  used  to  extensively  to  describe  such  domain,  where  shallow  water  and  its  depth,  as 
well  as,  shoreline  and  possible  fresh  water  composition  and  tidal  waves  play  a  critical 
role  in  battle  space  management  and  the  type  of  vessels  or  forces  that  USN  can  deploy  to 
this  domain. 


Figure  23.  Littoral  Ocean  Domain,  from  [16] 


The  term  “Littoral”  and  “Shallow”  water  domains  are  synonymous  for  this  study 
and  any  body  of  water  not  of  open-ocean,  less  than  or  equal  to  500  feet  depth  and  close  to 
shore  is  in  this  category.  Figure  24  shows  the  vast  majority  of  ocean  domain  as  open 
ocean  and  from  a  depth  perspective,  the  LOD  is  only  a  fraction  of  water  domain  above 
the  continental  shelf.  Flowever,  the  vast  majority  of  economic  transaction  and  military 
conflicts  are  conducted  in  LOD  and  this  number  continues  to  grow  [17].  Taking  it 
further,  the  LOD  will  be  bounded  generally  anywhere  but  open  ocean  and  within  the 
domains  of  depth  ranges  from  150  to  300  feet  or  approximately  10  to  20  times  the 
average  draft  of  USN’s  Littoral  Combat  Ship  (LCS,  12.8  feet  /  3.9  meters)  and  presence 
of  tidal  currents,  gradual  decent  of  ocean  bottom  shelf  to  open  ocean  depth,  presence  of 
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sedimentary  debris  throughout  water  and  numerous  bottom  sediment  types  to  inelude 
clay,  sand  and  rocks  are  some  of  many  characteristics  that  LOD  embodies. 


Figure  24.  World  Ocean  Domain,  from  [18] 


A.  LITTORAL  OCEAN  DOMAIN  RELATED  TO  UNDEX 

To  further  simplify  LOD  UNDEX  problems,  only  two  types  of  bottom  sediments 
will  be  considered  for  this  study — clay  and  sand — for  its  reflectivity  or  absorption  of 
energy  when  induced  to  forces  incurred  during  UNDEX.  The  distinct  differences  in 
behavior  of  cohesive  clay  bottom  type  and  non-cohesive  sand  will  not  be  explored  during 
this  study,  rather  its  reflectivity  or  absorption  of  energy.  Sand  in  general  will  absorb 
more  energy  from  UNDEX  due  to  its  increased  porosity  between  sand  particles.  Also, 
the  angle  of  decline  or  incline  in  shallow  water  areas,  like  beaches  or  shoreline  coming  in 
from  the  deep  abyss  will  be  assumed  to  be  small  or  negligible.  Therefore,  only  the  depth 
and  side  boundary  conditions  will  be  explored  for  UNDEX  phenomenon  in  this  research. 

There  are  four  primary  forces  that  can  be  induced  to  a  surface  vessel  operating  in 
an  EOD.  They  are  initial  detonating  shock  pressure  wave,  the  “hammer”  shockwave 
created  by  the  BC  formation  and  dissipation,  the  bottom  bounce  (BB)  or  laterally 
reflected  (ER)  pressure  waves  and  rising  cyclic  gas  bubbles  due  to  explosive  materials. 
While  all  four  forces  are  similar  to  that  of  what  is  experienced  in  an  open  ocean,  the 
effect  of  bottom  bounce,  LR  pressure  wave  and  additional  shockwave  due  to  cavitation  in 
EOD  is  what  sets  this  domain  apart.  The  rise  and  effect  of  gas  bubble  is  purely 
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determined  by  the  size  of  initial  bubble  radius;  therefore,  depending  on  the  eharge  size 
and  initial  bubble  radius,  the  subsequent  cyelie  bubbles  may  or  may  not  exist.  The  most 
distinctive  feature  that  differs  in  LOD  BC  formation  or  dissipation  is  the  follow  on 
disorderly,  additional  expansion  and  collapse  of  BC.  For  extremely  shallow  depth  below 
100ft,  the  bottom-reflected  shock  or  BB  arrives  before  cavitation  collapse  starts,  and  as  a 
result,  it  interacts  with  BC  in  a  complex  manner,  narrowing  its  region  and  hastening  its 
collapse  to  its  original  state  [11].  Additionally,  it  induces  additional  cavitation  in  the 
vicinity  of  the  bubble  and  gives  rise  to  a  more  disorderly  collapse  involving  the  original 
cavitation  region,  as  well  as,  additionally  induced  cavitated  areas.  The  shallower  depths 
also  play  a  role  in  formation  of  and  subsequent  cyclic  behavior  in  the  detonating  bubble. 
Again,  if  the  depth  of  LOD  is  smaller  than  the  radius  of  first  period  cyclic  wave  of 
detonation,  the  debris  within  the  detonating  sphere  will  spew  out  in  its  first  expansion 
creating  a  one  large  plume  rather  than  several  subsequent  plumes. 

B,  UNDEX  BENCHMARK  PROBLEM  FOR  LITTORAL  OCEAN  DOMAIN 

Using  DYSMAS,  a  benchmark  UNDEX  problem  for  this  study  is  created  and 
validated  (Figure  25).  This  fluid  or  Eulerian  domain  contains  three  distinct  layers  of 
Eulerian  materials  to  include  clay  bottom,  sea  water  and  a  layer  of  air.  The  origin  of  axis 
lies  at  water-air  interface  right  above  the  charge,  and  the  depth  and  lateral  distances  of 
domain  are  varied  in  order  to  meet  specific  objectives  for  each  run  cases.  Both  charge 
and  target  locations  are  prescribed  at  some  point  in  the  domain,  (X,  Y,  Z).  In  addition, 
HBX-1  will  be  used  as  the  charge  for  the  remainder  of  this  study,  due  to  carefully 
documented  characteristics  of  its  composition  and  UNDEX  behaviors.  HBX-1,  also 
known  as  High  Blast  Explosive  was  developed  during  WWII  as  desensitized 
modification  of  Torpex  explosives  and  its  composition  includes  38%  TNT,  40%  RDX, 
Aluminum  Wax  and  Eecithin  [19].  RDX,  also  known  as  Research  Department  Explosive 
is  an  explosive  nitroamine  widely  used  in  military  and  industrial  applications. 

The  setting  of  boundary  conditions  for  each  run  was  of  critical  importance;  since, 
they  determined  the  approximated  interactions  between  various  interfaces  of  UNDEX 
geometry  and  the  way  shockwaves  and  energy  will  travel  accordingly.  The  ocean  bottom 
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and  lateral  boundary  (B1  though  B4),  such  as  channel  walls,  were  assumed  to  be  a  rigid 
body  initially  to  simplify  the  modeling  requirements.  For  DYSMAS,  there  are  essentially 
three  boundary  settings:  free,  wall  and  0  to  99%  reflectivity.  The  condition  “free” 
includes  gravity  correction,  relevant  for  the  z-direction  only.  The  condition  “wall” 
specifies  a  reflecting  or  rigid  boundary  condition.  Lastly,  the  numbers  between  0  and  1 
can  be  used  for  partially  reflecting  boundary  conditions  [12]. 
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Figure  25.  UNDEX  Benchmark  Problem,  Defining  Littoral  Ocean  Domain 

To  confirm  initial  assumptions  made  regarding  LOD  in  reference  to  UNDEX 
parameters,  several  initial  case  studies  were  conducted  with  DYSMAS  by  first  simply 
varying  ocean  depths.  Initial  depths  were  varied  from  1000ft  to  100ft,  and  additional 
depth  ranges  were  further  broken  down  and  simulated  to  observe  or  validate  UNDEX 
characteristics  within  that  range.  The  simulated  cases  from  1  through  14  are  summarized 
below  in  Table  1.  The  UNDEX  characteristics  were  analyzed  using  several  results 
obtained  from  each  run  to  include:  flow  field,  shockwave  characteristics,  vertical  take¬ 
off  velocity  and  BC. 
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Case 

# 

Run 

Mode 

Euler 

Dimesions 

(X,Y,Z)  [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

1 

2D,  Euler 

800, 1, 1000 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

2 

2D,  Euler 

800,  1,  900 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

3 

2D,  Euler 

800,  1,  800 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

4 

2D,  Euler 

800,  1,  700 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

5 

2D,  Euler 

800, 1,  600 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

6 

2D,  Euler 

800, 1,  500 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

7 

2D,  Euler 

800,  1,  400 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

8 

2D,  Euler 

800, 1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

9 

2D,  Euler 

800,  1,  200 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

10 

2D,  Euler 

800, 1,  175 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

11 

2D,  Euler 

800,  1,  150 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

12 

2D,  Euler 

800, 1,  125 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

13 

2D,  Euler 

800,  1,  100 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

14 

2D,  Euler 

800, 1,  75 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Table  1.  Case  Studies  for  Defining  Littoral  Oeean  Domain 


1,  Flow  Field  Analysis 

The  overall  flow  field  results  for  deep  (greater  than  500ft)  and  shallow  (less  than 
300ft)  ocean,  or  LOD,  are  depicted  in  figures  26  and  27.  Both  flow  fields  depict  pressure 
contour  plots  or  “flow  field”  at  time  steps  of  (a)  10msec  (b)  35msec  (c)  50msec  and  (d) 
75msec.  These  time  steps  were  chosen  to  coincide  with  the  most  significant  UNDEX 
events  taking  place.  Both  figures  26  and  27  at  time  (a)  10msec  show  identical  initial 
shockwave  expansion  as  it  intercepts  the  water-air  interface.  Once  in  contact,  the 
characteristic  of  both  cases  take  a  different  turn.  The  expansion  and  initiation  BC  can  be 
clearly  seen  in  subsequent  field  plots  at  time  35  and  50msec,  with  initial  nucleation  of 
collapse  (Figure  26(b))  and  hammer  shockwave  (Figure  26(c)).  The  traveling  hammer 
shockwave  can  also  be  seen  in  Figure  26(d)  as  it  propagates  deep  into  the  ocean.  At  this 
point,  the  same  hammer  shockwave  have  already  reached  the  surface  and  dissipated.  As 
mentioned  in  previous  chapter,  hammer  shockwave  is  magnitude  of  the  pressure  pulse 
generated  by  the  collision  of  the  accreting  lower  and  upper  cavitation  boundary  at 
closure. 
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Figure  26.  500ft  Ocean  Depth  at  Time  Step  (a)  10  (b)  35  (c)  50  (d)  75msec 


The  flow  field  results  for  200ft  ocean  depth  are  depicted  in  Figure  27.  There  are 
several  contrasting  differences  to  previous  500ft  depth  flow  field  results.  First,  the 
bottom  bounced  (BB)  reflected  shockwave  are  shown  in  Figure  27(b).  The  magnitude  of 
this  BB  shockwave  at  reflection  is  identical  to  initially  radiating  shockwave.  As  it  travels 
back  up  towards  the  surface,  the  BB  shockwave  intercepts  the  “hammer”  shockwave 
from  collapsing  BC,  shown  in  Figure  27(c).  Last,  as  BB  shockwave  passes  through  the 
hammer  shockwave  and  hit  the  water-air  interface,  it  invokes  additional  BC  that  can  be 
seen  in  Figure  27(d).  Clearly  the  chaotic  nature  of  shallower  depth  UNDEX  can  be 
contrasted  from  the  500ft  flow  field  results  and  start  to  emerge. 
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Figure  27.  200ft  Ocean  Depth  at  Time  Step  (a)  10  (b)  35  (c)  50  (d)  75msec 


Again,  the  500ft  depth  UNDEX  results  were  representative  of  all  depth  cases 
deeper  than  300ft  and  200ft  depth  results  were  representative  of  cases  shallower  than 
300ft.  From  these  flow  field  results,  there  are  four  differing  characteristics  of  deep  and 
shallow  water  UNDEX.  First,  the  way  in  which  initial  shockwave  travels  and  interacts 
with  “other”  reflecting  waves  differ.  Unlike  open  ocean,  as  initial  shockwave  reflects  off 
shallow  bottom,  EB  or  water-air  interface  and  heads  towards  a  point  of  intersection  of  all 
three  waves,  additional  cavitation  and  shockwaves  are  created.  Second,  the  BC  collapse 
of  shallow  water  is  more  violent  or  higher  in  pressure  fluctuations  during  “hammer” 
shockwave  propagation  due  to  additionally  induced  cavitation  or  shockwaves.  Third,  the 
BB  shockwave  during  shallow  water  UNDEX  induces  additional  energy  back  into  the 
fluid  field  creating  additional  BCs  and  shockwaves  that  can  be  as  close  as  the  initial 
shockwave  and  energy.  The  additional  BC  can  be  seen  in  Figure  27(d)  as  the  BB 
shockwave  crashes  into  the  water-air  interface,  inducing  another  set  of  spallation  and 
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additional  formation  of  cavitations.  Lastly,  the  additional  BCs  that  are  formed  during 
BB’s  collision  of  water-air  interface  create  an  additional  “hammer”  shockwave  that  also 
can  be  as  high  as  the  initial  closure  pulse.  Another  flow  field  results  are  shown  in  Figure 
28  for  shallow  water  UNDEX  to  show  the  violent  and  multiple  cavitations  forming  due  to 
interactions  of  several  reflected  shockwaves  from  the  bottom.  As  shown,  the  white 
portion  of  flow  field  indicates  BC  formations.  Clearly,  the  additionally  induced  BC 
formation  and  its  chaotic  behavior  can  be  seen  and  at  times  greater  than  that  of  initial  BC 
in  size  and  scope. 


c.  +  51.58  ms  d.  +  121.6  ms 

Figure  28.  UNDEX  of  1 1021bs  TNT  at  Charge  Depth  of  66ft,  Water  Depth  of 

131ft,  from  [11] 


2,  Target  Pressure  Analysis 

The  open  ocean  and  EOD’s  target  pressure  profile  during  UNDEX  are  depieted  in 

figures  29  and  30.  Both  figures  show  that  the  initial  shoekwave  pressure  at  target 

location  to  be  approximately  55psi.  In  Figure  29,  this  initial  shoekwave  pressure  is 

followed  by  an  immediate  cutoff  pressure  at  8psi  with  reverberation  within  the  layer  of 

target  depth  and  air-water  interface,  and  finally  a  seeond  shoekwave  peak  pressure  at 

approximately  22psi  representing  the  “hammer”  shockwave  pressure  during  the  initial 

BC  closure.  To  get  a  perspective  on  shockwave  pressures,  14.7  psi  is  the  atmospheric 

pressure  that  is  felt  at  sea  level  and  55  psi  of  initial  peak  pressure  is  from  a  251bs  charge 

of  HBX-1  at  a  depth  of  50ft  of  sea  water.  55psi  is  elose  to  what  you  will  feel  at 
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approximately  120ft  deep  water,  stationary.  The  rest  of  the  time  period  show  negligible 
pressure  peaks  compared  to  initial  and  hammer  shockwave,  although  are  clearly  present. 
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Figure  29.  Target  Pressure  Profile  (300  to  1 000ft  Ocean  Depth) 


As  seen  on  the  overall  flow  field  plots  previously.  Figure  30  shows  additional 
pressure  waves  generated  due  to  BB  waves  and  additional  BC  closures.  For  75,  100  and 
125ft  of  water  depth,  the  additional  shockwave  present  due  to  BB  and  additional  BC  are 
found  to  be  up  to  90%  of  what  initial  shockwave  experienced.  For  the  100ft  case,  the  BC 
collapse  pulse  and  BB  SW  that  follow  are  almost  identical.  The  pressure  fluctuations  are 
clearly  more  observable  throughout  the  course  of  the  shallow  pressure  time  domain.  The 
pressure  fluctuations  felt  as  individuals  may  or  may  not  cause  structural  damage 
depending  on  the  magnitude  and  material  properties  of  structures  present.  However, 
accumulation  of  dynamic  nature  of  additional  shockwave  can  and  will  increase  the 
likelihood  of  the  structural  damage  due  to  dynamic  loads. 
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Target  Pressure  Profile  (75  ~  200ft  Ocean  Depth) 


Figure  30.  Target  Pressure  Profile  (75  to  200ft  Oeean  Depth) 

3,  Vertical  Take-Off  Velocity  Analysis 

Next,  the  vertieal  take-off  veloeities  (VTO)  are  analyzed  for  the  open  oeean  and 
shallow  water  UNDEX  events  (figures  31  and  32).  In  short,  VTO  represents  magnitude 
of  z-direction  veloeity  experieneed  at  target  location  during  UNDEX.  The  abrupt 
increase  in  velocity  at  a  constant  rate  followed  by  a  varying  decay  of  magnitude  shows 
close  resemblance  to  initial  shockwave  characteristics.  The  VTO  for  ocean  depths  of  300 
to  1000ft  is  shown  in  Figure  31  while  Figure  32  depicts  VTO  for  ocean  depth  of 
shallower  depths,  less  than  300ft.  Both  VTO  results  show  the  peak  of  velocity  at  1.75ft/s 
before  decaying  back  close  to  its  stationary  state  or  into  negative  velocity  due  to 
additional  fluid  movements  created  by  additionally  reflected  shockwaves  or  formed 
cavitations.  For  deeper  water,  rest  of  the  remaining  VTO  during  period  of  UNDEX 
seems  to  decrease,  approaching  close  to  negative  2ft/s. 
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Vertical  Take-Off  (VTO)  Velocity  (300  -  tOOOft  Ocean  Depth) 


Figure  3 1 .  Vertical  Take-Off  Velocity  (75  to  200ft  Ocean  Depth) 


An  entirely  different  phenomenon  takes  place  for  shallower  depth  in  Figure  32. 
Although  the  peak  VTO  is  identical  as  mentioned,  the  additional  VTOs  felt  due  to  BB 
and  additional  BC  formation  and  collapse  can  be  seen.  Again,  closely  resembling 
shockwave  pressure,  the  addition  VTO  experienced  at  target  locations  is  found  to  be  up 
to  and  even  as  greater  than  50%  of  the  initial  VTO.  The  numerous  additional  VTO  peaks 
and  chaotic  nature  of  shallower  depth  is  apparent  and  closely  resembles  previously 
analyzed  flow  fields  and  shockwaves.  In  short,  shallower  the  depth,  VTO  is  felt  more 
violently  and  duration  is  longer.  As  observed  in  pressure  profiles  previously,  the  distinct 
characteristics  of  less  than  300ft  starts  to  emerge  for  both  UNDEX  parameters. 
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Vertical  Take-Off  (VTO)  Velocity  (75  ~  200ft  Ocean  Depth) 


Figure  32.  Vertical  Take-Off  Velocity  (75  to  200ft  Ocean  Depth) 


4,  Bulk  Cavitation  Analysis 

Lastly,  the  BC’s  overall  characteristics  and  size  are  analyzed.  The  total  volume  of 
BC  formed  during  open  ocean  and  shallow  water  UNDEX  are  shown  in  figures  33 
through  35.  The  open  ocean  cases,  Figure  33,  show  a  peak  BC  volume  of  approximately 
6x10  ft  (4.5  X  10  Gallons).  This  is  a  tremendous  amount  of  water  moved  by  a  251bs 
charge  of  HBX-1;  hence,  like  the  initial  shockwave  of  55  psi,  one  can  gauge  at  the 
tremendous  energy  that  is  imparted  to  the  water  during  just  the  initial  part  of  an  UNDEX 
event.  For  all  300  to  1000ft  depth  oceans,  or  cases  without  bottom  reflection,  the  peak 
and  total  BC  volumes  are  almost  identical.  By  35msec,  BC  starts  to  collapse  and  by  50  to 
60msec,  complete  dissipation  of  initial  BC  can  be  observed  both  in  flow  field  plots 
previously  and  BC  volume  plots. 


40 


Bulk  Cavitation  Volume  (300  ~  1000ft  Ocean  Depth) 


Figure  33.  Bulk  Cavitation  Volume  (300  to  1000ft  Ocean  Depth) 


The  shallow  water  BC  characteristics  paint  a  different  story  as  briefly  mentioned 
in  previous  section.  The  increased  cavitated  zones  induced  by  the  BB  wave  and  are 
clearly  evident  in  shallow  water  (figures  34  and  35).  From  75  to  200ft  of  water,  the  total 
BC  formed  in  75ft  of  water  surpasses  all  other  depth’s  BC  formed  individually,  as  shown 
with  red  dotted  horizontal  line  in  Figure  34.  In  100ft  of  water,  the  BC  peak  is  observed  to 
be  6  X  10  ft  ,  over  four  million  gallons.  Four  million  gallons  is  approximately  equivalent 
size  of  present  day  water  tower  that  exists  throughout  various  cities.  This  dramatic 
increase  in  additional  BC  formation  for  shallower  depth  can  be  attributed  to  many  factors, 
but  for  this  study,  focus  will  be  on  characterizing  shallow  depth  and  effects  of  boundaries 
at  the  bottom  and  lateral  position  of  water  domains  that  are  present  in  shallow  ocean  or 
LOD. 
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X  10*  Cavitation  Volume  (75  ~  200ft  Ocean  Depth) 


Figure  34.  Bulk  Cavitation  Volume  (75  to  200ft  Ocean  Depth) 


X  10®  Bulk  Cavitation  Volume  (100  ~  200ft  Ocean  Depth) 


Figure  35.  Bulk  Cavitation  Volume  (100  to  200ft  Ocean  Depth  Cases) 
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V.  INVESTIGATION  OF  CHARGE  SIZE  AND  DEPTH 


Before  moving  further,  the  effects  of  charge  size  and  depth  must  be  addressed  and 
validated  for  LOD.  To  do  so,  the  same  UNDEX  fluid  domain  was  used  as  in  previous 
cases,  except  the  ocean  depth  was  maintained  at  500ft  deep  while  charge  sizes  and  depths 
were  varied.  The  depth  at  500ft  was  chosen  to  analyze  ocean  depth  of  greater  than  300ft, 
but  not  quite  deep  enough  to  see  the  negligible  changes  in  UNDEX  characteristics  caused 
by  charge  size  and  depth.  The  ocean  domain  for  this  case  studies  are  shown  in  Eigure 
36  and  Table  2,  where  cases  15  through  24  are  summarized  for  its  detailed  set  up. 


Z 


B4 


B3 


Eigure  36.  Littoral  Ocean  Domain  for  Charge  Size  and  Depth  Case  Studies 


Case 

# 

Run 

Mode 

Euler 

Dimesions 
(X,Y,Z)  [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

15 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

16 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

200 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

17 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

300 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

18 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

400 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

19 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

500 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

20 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

100 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

21 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

200 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

22 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

300 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

23 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

400 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

24 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

500 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Table  2.  Case  Studies  for  Varying  Charge  Size  and  Depth 


43 


A. 


VARYING  CHARGE  SIZE  (100  -  500LBS  OF  HBX-1) 


1,  Flow  Field  Analysis 

The  overall  flow  field  result  of  varying  charge  sizes  are  shown  in  Figure  37.  The 
time  step  chosen  represents  three  distinctly  important  events  taking  places:  (a) 
propagation  of  initial  shockwave  (b)  maximum  radius  of  BC  and  (c)  BC  collapse  where 
hammer  shockwave  is  generated.  Again,  as  the  initial  shockwave  slams  into  the  water-air 
interface,  the  BC  develops,  maturing  into  maximum  radius  by  55msec  and  eventually 
collapse  at  75msec  as  the  initial  shockwave  continues  to  propagate  into  the  depth  and 
remaining  ocean  domain.  The  chosen  500ft  ocean  depth  does  not  interfere  with  initiation 
and  collapse  of  BC  and  secondary  and  tertiary  UNDEX  events  are  not  readily  observed 
during  the  remainder  of  chosen  time  step  of  150msec,  hence  not  included.  For  rest  of  the 
charge  size  cases,  the  differences  in  the  flow  field  results  were  maximum  BC  radius  and 
time  of  BC  closure,  obviously  due  to  its  sizes. 


le.  ^  ^  1 

> 

(a) 

..M 

1 

L 

_ _ ^ 

Figure  37.  Flow  Field  at  Time  Step  (a)  10  (b)  55  (c)  75msec 


2,  Target  Pressure  Analysis 

The  results  of  target  pressure  for  varying  charge  sizes  are  depicted  in  figures  38 
and  39.  Figure  38  shows  the  maximum  pressure  of  5001bs  FlBX-1  at  approximately 
l,000psi.  Next,  Figure  39  shows  maximum  pressure  peaks  of  100  to  4001bs  HBX-1.  The 
difference  of  peak  pressures  between  400  and  5001bs  HBX-1  is  close  to  700psi,  while 
differences  of  peak  pressure  between  100  and  4001bs  HBX-1  is  close  to  200psi.  These 
drastic  and  almost  exponential  differences  in  pressure  peak  between  500  and  4001bs 
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charge  is  of  significant,  since  it  ean  be  used  to  bound  problems  via  similitude  to  pick  an 
ideal  weight  that  is  most  suitable  for  a  given  case  studies. 


Target  Pressure  for  Various  Charge  Sizes 


Figure  38.  Target  Pressure  Profile  (100  to  SOOlbs  Charge  Weight) 


Target  Pressure  for  Various  Charge  Sizes 


Figure  39.  Initial  Target  Pressure  Profile  (100  to  SOOlbs  Charge  Weight) 
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3,  Vertical  Take-Off  Velocity  Analysis 

The  VTO  results  for  various  eharge  sizes  are  depicted  in  Figure  40.  The 
maximum  VTO  for  500Ibs  peaks  at  approximately  37ft/s  while  4001bs  or  less  charges 
peak  at  12ft/s  or  less.  Similar  to  pressure  peaks  observed  previously,  the  differences  in 
the  VTO  between  100  to  4001bs  charges  and  SOOlbs  charges  are  quite  drastic.  Other 
characteristics  that  are  noticeable  are  minor  VTO  that  starts  to  take  shape  after  50msec 
for  4001bs  or  less  charge  sizes.  These  minor  VTO  can  be  attributed  to  the  smaller  BC 
collapsing  and  inducing  additional  fluid  movements  that  can  be  observed  during  the  time 
step  of  simulations.  For  5001bs  charge,  the  sheer  size  fluids  moved  and  BC  formed 
during  detonation  takes  a  lot  longer  to  return  to  its  original  state  than  the  designated 
150msec  time  frame,  hence,  the  magnitude  of  its  VTO. 


Vertical  Take-Off  Velocity  (VTO)  for  Various  Charge  Sizes 


Figure  40.  Vertical  Take-Off  Velocity  (100  to  5001bs  Charge  Weight) 


4,  Bulk  Cavitation  Analysis 

The  BC  formation  and  dissipation  for  various  charge  weight  is  shown  in  Figure 
41.  Again,  the  drastic  differences  in  BC  volume  between  100  to  4001bs  charge  weight 
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compared  to  SOOlbs  charge  cannot  be  ignored.  The  shape  and  speed  of  cavitations  are 
consistent  with  previous  cases  studied.  The  sheer  amount  of  water  that  is  cavitated  due  to 
energy  imparted  into  the  water  via  propagating  initial  shockwave  is  close  to  3x10  ft 
(2x10^  Gallons)  for  5001bs  charge,  while  4001bs  charge  cavitates  close  to  7.5x10^  ft^ 
(6x10^  Gallons)  water  domain.  It  is  important  to  keep  in  mind  that  at  this  point,  the 
effect  of  BB  or  LB’s  reflected  waves  and  their  cause  of  additional  cavitation  are  not  even 
considered  due  to  chosen  500ft  ocean  depth.  At  this  depth,  for  the  chosen  charge  weight, 
the  BC  collapses  prior  to  any  interfering  additional  reflected  shockwave  present. 


X  10^  Caviation  Volume  for  Various  Charge  Sizes 
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Figure  41 .  Bulk  Cavitation  Volume  (100  to  5001bs  Charge  Weight) 


B.  VARYING  CHARGE  DEPTH  (100  -  500FT  OF  lOOLBS  HBX-1) 

1,  Flow  Field  Analysis 

The  overall  flow  field  result  of  varying  charge  depths  are  shown  in  figures  42  and 
43.  Figure  42  depicts  a  300ft  charge  depth  case  and  Figure  43  depicts  500ft  charge  depth 
case  where  the  charge  itself  is  at  the  water-clay  (bottom)  interface.  The  time  step  chosen 
to  represent  four  distinctly  important  events  taking  places  are;  (a)  propagation  of  initial 
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shockwave  (b)  initial  shockwave  reaching  water-air  interface  (c)  maximum  BC  radius 
and  (d)  additional  cavitation  formed  due  to  BB  shockwaves. 


^  (c)  p""  (d) 

Figure  42.  Flow  field  at  Time  Step  (a)  10  (b)  60  (c)  73  (d)  150msec  (300ft 

Charge  Depth) 

Unlike  previous  cases,  the  shockwave  in  Figure  42  takes  a  lot  longer  to  reach  the 
water-air  interface  due  to  its  charge  depth.  As  a  result,  the  overall  pressure,  VTO  and  BC 
events  are  delayed  up  to  60msec.  Further,  the  shockwave  pressure  and  its  magnitude  that 
reached  the  water-air  interface  has  drastically  reduced,  shown  in  Figure  42(b),  since  its 
energy  has  been  continuously  dissipating  into  the  water  as  it  travels  up  to  the  surface. 
The  usual  BB  shockwaves  that  are  reflected  passes  the  charge  depth  (Figure  42(c))  and 
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reaches  water-air  interface  after  the  BC  have  collapsed.  The  collapse  of  additional 
cavitation  formed  by  BB  shockwave  can  be  seen  in  Figure  42(d).  Next,  for  500ft  charge 
depth,  two  distinctly  different  events  take  place,  shown  in  Figure  43.  First,  since  the 
charge  is  embedded  into  the  bottom  clay,  as  detonation  takes  place,  the  initial  shockwave 
and  BB  shockwaves  are  merged  to  form  a  single  shockwave  that  propagates  up  to  the 
surface  of  the  water.  As  this  emerging  shockwave  travels  up  to  the  water-air  interface 
and  collides,  the  BC  formed  is  much  greater  than  any  other  charge  depth  combined. 


Figure  43.  Flow  field  at  Time  Step  (a)  10  (b)  102  (c)  119  and  (d)  130msec 

(500ft  Charge  Depth) 


49 


From  the  flow  field  results,  this  huge  BC  ean  be  seen  in  Figure  43(e)  and  (d).  As 
this  BC  eontinues  to  grow,  a  similar  oharaeteristie  ean  be  observed  in  terms  of  additional 
eavitation  formed  during  shallow  water  eases  observed  previously.  In  short,  due  to  a 
eombined  shoekwave  in  lieu  of  two  initial  and  BB  shoekwaves,  the  BC  formed  for  this 
ease  is  almost  double  in  its  original  size  if  the  eharge  depth  was  deeper  than  its  half-way 
point  of  water  domain  or  submerged  into  the  bottom  elay.  For  eharge  depth  eloser  to 
water-air  interface,  similar  phenomenon  can  be  observed,  except  the  absence  of  BC  due 
to  50%  of  detonation  energy  expanded  directly  into  the  atmosphere. 


2,  Target  Pressure  Analysis 


Target  Pressure  for  Various  Charge  Depth 


Figure  44.  Target  Pressure  Profile  (100  to  500ft  Charge  Depth) 


The  overall  target  pressure  profiles  are  depicted  in  Figure  44  for  various  charge 
depth  cases.  The  pressure  peak  of  100ft  charge  depth  shows  maximum  pressure  of 
approximately  135psi  and  500ft  charge  depth  shows  a  maximum  pressure  peak  at 
approximately  40psi.  As  mentioned  previously  in  overall  flow  field  results,  as  the  charge 
depth  is  increased  and  with  it,  the  distance  to  water-air  interface  is  also  increased;  hence, 
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the  actual  magnitude  of  shockwave’s  energy  that  collides  with  the  interface  is  greatly 
reduced  since  it  must  travel  through  greater  volume  of  water.  As  a  result,  the  total 
volume  and  timing  of  BC  formation  is  greatly  reduced  and  delayed,  except  for  500ft 
charge  depth  case.  Taking  closer  look  at  the  500ft  charge  depth  case,  the  secondary 
shockwave  is  much  greater  than  that  of  its  counterpart  (Figure  45).  As  described 
previously  through  overall  flow  field  results,  this  is  due  to  the  emergence  of  two 
shockwaves:  initial  shockwave  and  BB  shockwaves.  This  secondary  shockwave  is 
similar  in  magnitude  of  lOOft  charge  depth  cases.  The  tertiary  and  rest  of  the  shockwaves 
observed  are  consistent  in  magnitude  and  time  duration  except  for  500ft  charge  depth 
case.  For  this  case,  tertiary  and  rest  of  the  shockwaves  are  actually  missing  since  these 
shockwaves  were  combined  to  form  greater  secondary  shockwave. 


Target  Pressure  for  Various  Charge  Depth 
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Figure  45.  Close-Up  Target  Pressure  Profile  (100  to  500ft  Charge  Depth) 


3,  Vertical  Take-Off  Velocity  Analysis 

The  VTO  for  various  charge  depth  cases  are  depicted  in  Figure  46.  Consistent 
with  previous  UNDEX  parameters,  the  decrease  in  VTO  as  charge  depth  is  increased  can 
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be  observed.  The  100ft  charge  depth’s  VTO  is  close  to  5ft/s  while  500ft  charge  depth  is 
shown  to  exhibit  VTO  of  close  to  1.75ft/s.  Not  surprisingly,  the  differences  to  the  400ft 
charge  depth’s  VTO  is  minimal  when  compared  back  to  the  bottom  depth  case  of  500ft. 
Another  characteristic  of  VTO  that  can  be  observed  is  that  as  the  initial  VTO  comes  back 
to  the  original  state,  this  overall  state  is  increased  for  shallower  depth.  For  example, 
VTO  for  100ft  returns  back  to  greater  than  0.5ft/s  and  hovers  at  this  state  of  velocity 
through  the  remainder  of  its  timeframe.  This  characteristic  is  also  true  for  the  remaining 
cases,  except  for  the  bottom  depth  case,  as  pointed  out  in  Figure  44.  Instead  of  initial 
VTO  spike  and  returning  to  constant  secondary  state,  there  is  a  secondary  spike  in  VTO 
followed  by  tertiary  fluctuated  spikes.  Again,  these  secondary  and  tertiary  spikes  are 
attributed  to  the  additional  cavitation  formation  at  initial  shockwave  interception  of 
water-air  interface,  rather  than  later  time  for  rest  of  the  depth  cases.  Lastly,  the  tertiary 
VTO  spike  settles  out  at  a  lower  magnitude  than  what  the  previous  cases  exhibited. 


Vertical  Take-Off  Velocity  (VTO)  for  Various  Charge  Depth 
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Figure  46.  Vertical  Take-Off  Velocity  (100  to  500ft  Charge  Depth) 
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4,  Bulk  Cavitation  Analysis 

The  BC  volume  formation  is  depicted  in  Figure  47.  The  overall  formation  and 
dissipation  of  BC  is  consistent  with  what  was  observed  in  previous  cases  where  BC  is 
initiated  at  initial  shockwave’s  interception  of  water- air  interface,  reaches  its  maximum 
radius  at  some  point  depending  on  the  charge  size  and  additional  cavitations  are  formed 
depending  on  the  presence  or  lack  of  bottom.  The  BC  volume  for  charge  depth  of  100 
and  200ft  follows  this  consistent  pattern  observed  previously;  however,  the  rest  of  the 
cases  differ.  For  a  charge  depth  of  300  and  400ft,  an  additional  “hump”  can  be  observed 
at  the  posterior  section  of  BC  formation  and  dissipation.  This  is  due  to  BB  shockwave 
catching  up  to  the  BC  dissipation  and  as  it  interferes  with  the  BC  collapse,  the  BB 
shockwave  induces  additional  cavitation  to  be  formed.  This  phenomenon  is  greatly 
amplified  in  the  500ft  charge  depth  case  where  initial  BC  and  additional  cavitation  are 
formed  simultaneously,  causing  greater  volume  of  overall  BC.  As  observed,  the  overall 
max  BC  formed  for  the  bottom  depth  case  is  close  to  2.3x10  ft  (1.7x10  Gallons),  while 
the  100ft  charge  depth  case’s  overall  max  BC  volume  is  drastically  less  than  the  former  at 
1.8x10^  ft^  (1.3x10^  Gallons).  While  the  amount  of  overall  BC  formed  is  not  as 
important  as  additional  cavitations  formed  due  to  additional  shockwaves  that  can  be 
created,  the  sheer  volume  of  additional  BC  created  due  to  charge  depth  differences  cannot 
be  ignored.  Hence,  the  collapse  of  BC  formed  by  the  500ft  charge  depth  will  have 
greater  residual  shockwaves  originated  from  its  hammer  pulse. 
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X  10^  Caviation  Volume  for  Various  Charge  Depth 


Figure  47.  Bulk  Cavitation  Volume  (100  to  500ft  Charge  Depth) 
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VI.  INVESTIGATION  OF  LATERAL  BOUNDARY  CONDITIONS 


In  LOD,  the  lateral  boundary  (LB)  eondition  (B2)  play  a  role  in  truncating  BC 
formation,  as  well  as,  drastically  enhancing  the  role  of  reflected  wave,  much  like  the  BB 
effect  observed  previously.  Cases  such  as  channels,  test  ponds  and  areas  of  LOD  where 
lateral  spaces  are  limited,  are  critical  to  in-depth  analysis  of  this  case  due  to  its  effects  on 
UNDEX  parameters.  To  take  a  look  at  such  phenomenon,  the  placement  and  reflectivity 
of  lateral  wall  boundaries  were  varied  during  UNDEX  (Eigure  48).  The  placement  of  EB 
conditions  (B2.1,  2.2,  and  2.3)  were  a  matter  of  where  the  BC’s  edge  was  located  and 
accordingly  placed  inside,  at  the  edge  or  outside  of  the  maximum  BC  radius.  These 
represent  the  physical  boundaries  that  the  target  and  charge  placements  are  as  prescribed 
in  the  Case  25  to  28  in  Table  3.  In  addition,  reflectivity  of  boundary  case  26  (B2.1), 
within  the  maximum  BC  radius,  was  varied  to  analyze  its  effect  and  influences  on 
UNDEX  parameters. 


z 


Eigure  48.  Eateral  Wall  Boundary  Conditions  (B2. 1,  B2.2  and  B2.3) 
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Case 

# 

Run 

Mode 

Euler 

Dimesions 

(X,Y,2)  [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

25 

2D,  Euler 

230,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

26 

2D,  Euler 

164,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

26.1 

2D,  Euler 

164,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

0.50 

Wall 

Wall 

0.656 

26.2 

2D,  Euler 

164,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

27 

2D,  Euler 

197,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

28 

2D,  Euler 

230,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

Table  3.  Case  Studies  for  Varying  Lateral  Boundary 


A,  VARYING  LOCATIONS 

The  next  set  of  numerical  results  show  initial  and  follow  on  shockwave  pressure, 
vertical  take-off  velocity  and  BC  volume  experienced  at  target  location  and  depth.  As 
mentioned  in  Table  3  for  case  studies  25  through  28,  the  lateral  target  location  is  at  100ft 
and  0.820ft  deep,  for  all  practical  purpose  in  regards  to  mesh  size  at  0.656ft,  this  depth  is 
at  shallow  surface  depth  near  water-air  interface. 

1,  Flow  Field  Analysis 

The  overall  pressure  flow  field  results  are  shown  in  figures  49  to  52.  Figure  49  is 
the  base  case  where  LB  was  set  well  outside  the  maximum  BC  radius  with  free  boundary 
condition.  The  snapshots  of  flow  field  were  taken  at  time  step  10,  30,  50  and  111msec 
where  significant  UNDEX  events  took  place.  Consistent  with  previous  cases,  the  flow 
field  results  show  expected  behaviors  of  (a)  radiating  initial  shockwave,  (b)  initiation  and 
dissipation  of  BC,  (c)  BC’s  “hammer”  shockwaves  and  (d)  influence  of  other  reflected 
waves  from  bottom  or  lateral  boundaries.  In  Figure  49(d),  the  BB  SWs  are  shown 
passing  the  charge  depth  as  it  intercepts  the  water-air  interfaces.  The  “free”  LB  condition 
set  at  x-maximum  of  this  domain  seems  to  have  no  observable  effect  on  the  radiating 
shockwave  within  flow  field. 
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Figure  49.  Case  25  (Base),  Time  Step  (a)  10  (b)  30  (c)  50  (d)  1 1 1msec 

Next,  case  26  involved  setting  the  LB  conditions  inside  the  maximum  BC  radius 
to  observe  its  effect  on  UNDEX  during  the  identical  time  steps  of  case  25.  No 
differences  were  observed  at  time  step  10msec  where  initial  shockwave  intercepts  the 
water-air  interface,  shown  in  Figure  50(a).  However,  subsequent  time  steps  clearly  show 
distinction  on  the  effect  of  wall  boundary  condition  set  inside  the  BC  radius.  First, 
Figure  50(b)  shows  BC  at  a  maximum  radius  similar  to  that  of  the  base  case,  but  with 
shallower  maximum  depth.  Second,  Figure  50(c)  shows  the  presence  of  secondary 
shockwave  that  is  reflected  off  the  FB  as  it  travels  back  towards  the  charge  location, 
creating  turbulence  and  additional  cavitation  following  its  path  as  it  intercepts  the 
hammer  pressure  shockwaves.  At  this  point,  the  additional  cavitation  following  this 
wave  is  also  collapsing,  creating  another  shockwave,  although  smaller  in  magnitude.  As 
this  reflecting  wave  intercepts  the  left  boundary  condition  and  again  reflected  back  into 
the  water  domain,  the  BB  waves  are  seen  disrupting  this  process  as  it  begins  to  dissipate 
completely  post  time  step  of  1 1 1msec,  shown  in  Figure  50(d). 
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Figure  50.  Case  26  (Inside  BC),  Time  Step  (a)  10  (b)  30  (e)  50  (d)  1 1  Imsee 


Case  27  shows  the  LB  plaeement  at  the  edge  of  the  BC  radius  (Figure  51).  The 
flow  field  at  time  step  10  and  30msee  are  identieal  to  base  ease  in  both  shoekwave 
propagation  and  BC’s  max  radius  and  depth.  At  time  step  50msee  of  Figure  51(e),  the 
refieeted  LB  SW  begins  to  emerge  with  trails  of  additional  cavitation  forming.  This  LB 
shockwave  is  delayed  due  to  the  boundary’s  placement  and  the  magnitude  seems  to  be 
similar  to  that  of  case  26.  Lastly,  Figure  51(d)  shows  the  several  BB  and  LB  shockwaves 
intercepting  to  initiate  additional  turbulence  and  cavitation  formations. 
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Figure  51.  Case  27  (Edge  of  BC),  Time  Step  (a)  10  (b)  30  (c)  50  (d)  1 1 1msec 

The  results  of  LB  placement  well  outside  of  the  BC  radius  (case  28)  are  shown  in 
Figure  52.  The  flow  field  at  time  step  10  and  35msec  are  identical  to  the  base  case  in 
both  shockwave  propagation  and  max  radius  and  depth  of  BC.  At  time  step  50msec, 
shown  in  Figure  52(c),  the  delayed  LB  reflected  wave  have  yet  to  intercept  the  hammer 
pressure  wave  generated  as  initial  BC  collapses.  Once  again,  following  closely  to  the 
laterally  reflected  wave,  an  initiation  of  additional  cavitation  is  clearly  visible  in  Figure 
52(c).  Lastly,  intercepting  laterally  reflected  and  the  BB  SWs  are  shown  in  Figure  52(d). 
As  expected  this  interception  takes  place  at  a  time  step  longer  than  any  other  cases 
discussed  due  to  the  location  of  LB  conditions  and  compared  to  previous  cases  the 
additionally  induced  cavitation  seems  to  have  dissipated  completely  by  this  time  step  of 
111msec.  To  further  analyze  various  UNDEX  parameters  numerically,  a  point  target 
location  was  chosen  at  100ft  lateral  distance  and  surface  depth  from  the  charge  location. 
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2,  Target  Pressure  Analysis 

The  resulting  shockwaves  at  the  target  location  are  analyzed  for  various  LB 
conditions  (figures  53  to  55).  The  first  distinct  spikes  at  approximately  25msec  are  initial 
compressive  shocks  as  it  reaches  the  target,  immediately  followed  by  tensile  SW  as 
water-air  interface  expands  and  spallates  to  initiate  BC  formation,  the  BC  closure  or 
“hammer”  SW  and  followed  by  laterally  returning  reflected  shockwave’s  interception  of 
target  location.  The  BC  collapsing  shockwave  is  considered  to  be  secondary  and  laterally 
reflected  wave  is  considered  tertiary  during  this  case  study.  The  other  shockwaves 
observed  post  110msec  is  BB  SWs.  The  location  and  magnitude  of  these  SWs  are 
consistent  theoretically.  However,  the  timing  and  where  both  LB  and  BB  reflected  SW  is 
located  during  the  entire  pressure  time  step  can  cause  undue  damage  to  the  ship  by 
invoking  modal  response.  Another  characteristic  to  point  is  the  fact  that  the  magnitude  of 
the  LB  and  BB  shockwaves  are  at  time  even  greater  than  that  of  BC  collapsing  SWs. 
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Target  Pressure  for  Various  Lateral  Boundary  Conditions 


Figure  53.  Target  Pressure  for  Various  Lateral  Boundaries 


Target  Pressure  for  Various  Lateral  Boundary  Conditions 


Figure  54.  Initial,  Secondary  and  Tertiary  Target  Pressure  for  Various  Lateral 

Boundaries 
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The  initial  shockwave  magnitude  for  all  LB  cases  are  identical  in  magnitude  at 
approximately  55psi  and  shaped  similarly;  however,  what  follows  after  the  BC  collapse, 
differs  depending  on  the  size  of  the  water  domain  or  location  of  LB  conditions.  Between 
approximately  25  and  35msec  where  BC  collapses,  the  magnitude  of  cyclic  trapped 
shockwave  between  the  surface  and  target  location  is  largest  for  LB  condition  located 
outside  BC  max  radius  (Figure  54).  Next,  the  magnitude  of  BC  collapse  shockwave  also 
differs  with  this  boundary  condition  at  22  psi.  What  follow  next  are  three  tertiary 
pressure  peaks  showing  laterally  reflected  waves  for  all  three  boundary  conditions  with 
closest  boundary  showing  highest  and  earliest  magnitude  of  reflected  shockwave.  Lastly, 
the  group  of  tensile  and  compressive  shockwave  after  1 10msec  show  various  shockwaves 
as  both  BB  and  LB  reflected  shockwave  intersect  to  cause  additional  cavitations  and 
shockwaves  where  it  reaches  peak  magnitudes  of  greater  than  hammer  shockwave 
(Figure  55). 


Target  Pressure  for  Various  Lateral  Boundary  Conditions 


Figure  55.  Bottom  Bounce  Target  Pressure  for  Various  Boundaries 
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3,  V ertical  T ake-Off  V elocity  Analysis 

The  VTO  for  various  LB  conditions  are  shown  in  Figure  56.  The  initial  upward 
velocity  peaks  at  1.75ft/s  for  all  LB  conditions.  What  follows  next  mirror  the  pressure 
profile  seen  in  previously.  As  the  initial  peak  VTO  settles  to  near  zero  velocity,  the 
hammer  shockwave  at  approximately  at  37msec  induces  another  VTO  back  to  .25ft/s.  As 
VTO  reattempts  to  settle,  the  three  additional  peak  pressure  waves  induced  by  LB 
between  49  to  75msec  causes  three  additional  peaks  of  VTO.  The  VTO  settles  once 
again,  and  then  picks  up  approximately  120msec  when  collision  of  LB  and  BB  reflected 
SWs  cause  additional  BC  and  dynamic  motion  of  water  within  the  domain.  It  is 
important  to  note  that  VTO  does  not  return  below  zero  from  approximately  35  to 
130msec  and  continues  to  be  induced  with  additional  forces  and  its  velocity  caused  by 
BC,  collisions  of  LB  and  various  reflected  wave. 


Vertical  Take-Off  Velocity  (VTO)  for  Various  Lateral  Boundary  Conditions 


Figure  56.  VTO  for  Various  Lateral  Boundaries 
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4,  Bulk  Cavitation  Analysis 

The  BC  volumes  at  various  LB  eonditions  eompared  to  the  base  eondition  are 
depleted  in  Figure  57.  Observing  the  initial  BC  formation  and  dissipation  from  time 
steps  0  to  50msec,  there  are  several  distinct  differences  in  various  LB  conditions.  First, 
base  and  LB2.3  conditions  are  practically  identical  in  nature  except  there  is  a  distinct 
“hump”  (1.25x10  ft  )  initiated  at  50msec  and  also  post  135msec.  This  phenomenon,  also 
observed  in  Figure  52(c),  is  the  initiation  of  additional  cavitation  at  LB  at  time  of  initial 
shockwave  reflection.  This  “hump”  in  BC  volume  is  not  observed  in  base  case  obviously 
due  to  different  LB  condition  setting  as  “free”  condition.  Second,  looking  at  the  LB2.1 
and  LB2.2,  a  similar  hump  is  visible  at  an  earlier  time  step  (near  37msec)  due  to  closer  in 
proximity  of  the  wall.  Lastly,  the  overall  BC  volume  of  LB2.1  is  much  greater  than  rest 
compared  to  the  base  and  other  LB  cases.  Like  shallower  depth  BB  shockwave  inducing 
enormous  volume  of  additional  cavitation,  the  same  effects  are  observed  here  for  LB 
conditions  well  within  the  BC’s  max  radius. 


X  -10^  Caviation  Volume  for  Various  Lateral  Boundary  Conditions 


Figure  57.  Bulk  Cavitation  Volume  for  Various  Lateral  Boundaries 
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The  next  range  of  time  step  from  50  to  lOOmsee  shows  the  effect  of  LB’s 
interferences  and  additional  cavitation  formed  as  the  lateral  shockwave  intercepts  the 
hammer  shockwave.  Clearly,  the  magnitude  and  speed  of  the  additional  cavitation 
formed  are  much  higher  for  LB2.1.  The  post  1 00msec  ranges  show  drastic  increase  in 
additional  cavitation  formation  due  to  BB  inclusion  of  already  chaotic  water  domain. 
Consistent  in  behavior  to  the  50  to  1 00msec  range,  the  LB  that  is  the  closest  to  the  charge 
had  the  largest  additional  cavitation  formed  surpassing  the  initial  cavitation  volume  for 
all  cases.  It  is  surprising  to  note  that  the  LB  conditions  have  almost  the  same  effect  in 
inducing  additional  cavitation  as  the  BB  has  on  BC. 

B,  VARYING  REFLECTIVITY  INSIDE  THE  BC  RADIUS 

In  order  to  further  investigate  the  degree  of  shockwave  propagation  as  it  pertains 
to  the  percentage  of  reflectivity  of  the  LB  conditions,  case  26  (LB2.1)  were  simulated 
once  again  using  50%  and  free  boundary  conditions.  For  case  27  (LB2.2)  and  28  (LB2.3) 
the  effects  of  such  change  in  reflectivity  due  to  greater  distance,  in  reference  to  max  BC 
radius,  were  not  as  drastic  as  case  26.  All  other  parameters  were  maintained  as 
previously  described  in  Table  3. 

1.  Flow  Field  Analysis 

The  flow  field  results  for  50%  reflective  LBs  are  shown  in  Figure  58. 
Comparable  to  base  case  26,  where  boundary  condition  was  set  as  100%  reflective  (wall), 
the  general  features  of  radiating  initial  shockwave,  BC  formation,  reflected  shockwave 
and  its  interception  of  BB  waves  and  BC  collapse  wave  are  nearly  identical.  Case  26  is 
the  base  case  in  this  section  where  its  LB  conditions  is  set  as  wall  and  case  26.1  is  set  at 
50%  reflectivity.  As  expected,  the  magnitude  of  shockwave  for  the  50%  reflective  wall  is 
clearly  lower  than  the  case  26  but  its  shockwave  travels  at  a  pattern  consistent  with  case 
26.  The  flow  field  at  111msec  is  drastically  different  than  that  of  base  case  26.  As 
Figure  58(d)  shows,  the  magnitude  in  which  BB  wave  interacts  with  laterally  reflected 
waves  are  a  lot  higher  than  the  base  case.  Since,  50%  of  laterally  reflected  energy  was 
absorbed  by  the  LB;  the  intercepting  wave  has  less  energy  to  neutralize  BB  waves  as  it 
traveling  up  to  the  water-air  interface.  In  short,  as  BB  waves  travels  up  and  passes  the 
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charge  depth,  the  energy  dissipation  is  a  lot  lower  for  50%  lateral  wall  reflective  case 
26.1.  The  result  is  increased  UNDEX  forees  applied  to  object  in  its  path  and  on  surface 
of  water. 


Figure  58.  Case  26.1  (Inside  of  BC,  50%  Reflectivity),  Time  Steps  (a)  10  (b)  30 

(c)  50  (d)  111msec 

Next,  the  flow  field  result  for  Case  26.2  is  shown  in  Figure  59.  Case  26.2  depicts 
case  26  with  free  or  no  reflective  FB  conditions  well  within  the  max  radius  of  BC.  The 
flow  field  result  is  almost  identical  to  the  base  ease  25  for  all  chosen  time  step,  except  for 
111msec.  At  approximately  100ft  deep  below  charge  depth  there  is  a  pressure  spike 
where  BB  SWs  meet  the  downward  traveling  BC  hammer  SW  (Figure  59(d)).  As  seen, 
where  these  two  pressure  waves  collide,  there  is  a  noticeable  spike  that  is  similar  to  when 
the  laterally  refleeted  wave  and  BB  wave  intercept.  Although  this  seems  minor,  the 
distinct  and  unexpected  difference  in  pressure  spike  can  cause  additional  forces  applied  to 
the  floating  object  within  the  timeframe  of  UNDEX  event. 
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Figure  59.  Case  26.2  (Inside  of  BC,  No  Reflectivity),  Time  Steps  (a)  10  (b)  30 

(c)  50  (d)  111msec 

2,  Target  Pressure  Analysis 

The  pressure  profde  for  cases  25,  26,  26.1  and  26.2  are  shown  in  figures  60  to  62. 
Observing  these  pressure  profiles  of  four  different  cases,  the  tertiary  LB  SW  following 
BC  hammer  shockwave  is  clearly  absent  except  for  case  26  and  26.1  where  LBs  were  set 
as  a  rigid  wall  and  50%  reflective  wall.  The  observed  shockwaves  are  expected  for  both 
cases  in  terms  of  their  location  and  magnitude.  A  closer  look  at  timeframe  from  20  to 
50msec  (Figure  61)  shows  base  case  25  and  case  26.2  pressure  profiles  are  indeed 
identical  for  both  BC  collapse  shockwaves  of  17.5psi  and  absence  of  LB  reflection  spike 
near  48msec.  Likewise  for  case  26  and  26.1,  the  pressure  profiles  also  closely  resemble 
each  other,  being  identical  in  BC  collapse  pressure  spike  of  22psi  at  37msec  and  presence 
of  laterally  reflected  wave  near  48mse,  with  wall  LB  setting  having  the  higher  pressure 
spike  at  23psi. 
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Target  Pressure  for  Case  26  at  Various  Reflectivity 


Figure  60.  Target  Pressure  for  Varying  Boundary  Refleetivity 


Target  Pressure  for  Case  26  at  Various  Relectivity 


Figure  61 .  Initial  and  Secondary  Target  Pressure  for  Varying  Boundary 

Reflectivity 


68 


Taking  a  closer  look  at  time  step  ranges  from  100  to  ISOmsee  where  BB  SWs 
start  to  emerge,  again  the  presenee  or  lack  of  LB  conditions  play  a  crucial  role  in 
invoking  additional  shoekwaves  during  this  time  frame.  Both  cases  25  and  26.2  show 
similar  pressure  fluetuations  between  13.5  and  IVpsi.  Cases  26  and  26.1  show  pressure 
fluetuations  between  6  and  24psi.  The  differences  between  wall  and  free  LB  eonditions 
are  distinct  with  possibility  of  LB  condition  invoking  higher  pressure  fluetuations  and 
within  the  ehanges,  the  refleetivity  setting  closely  resembles  how  mueh  fluetuations  are 
observed  for  oases  26  and  26. 1 . 


Target  Pressure  for  Case  26  at  Various  Reflectivity 


Figure  62.  Bottom  Bounce  Target  Pressure  for  Varying  Boundary  Reflectivity 

3,  Vertical  Take-Off  Velocity  Analysis 

The  VTO  results  are  shown  in  Figures  63.  The  initial  spike  in  VTO  peaks  at 
L75ft/s  for  all  oases.  The  seoondary  or  LB  induoed  VTO  for  all  oases  are  different  to 
some  degree,  but  the  presenee  of  spike  is  apparent  for  presenee  of  LB  setting  at  50%  or 
wall.  For  oases  25  and  26.2,  the  spike  is  absent;  however,  there  is  a  minor  differenoe 
between  the  two.  Again  for  BB  induced  VTO,  there  are  distinot  differenoes  between 
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presence  of  LB  conditions  and  higher  magnitudes  are  observed  for  cases  where  50%  or 
wall  setting  is  present.  Lastly,  compared  to  cases  27  and  28  (Figure  56),  the  absence  of 
varying  LB  induced  spikes,  as  well  as,  reduced  fluctuation  in  BB  induced  VTO  are  also 
apparent  for  this  case  (Figure  63). 


Vertical  Take-Off  Velocity  (VTO)  for  Case  26  at  Various  ReftectMty 


Figure  63.  VTO  for  Various  Lateral  Boundary  Reflectivity 

4,  Bulk  Cavitation  Analysis 

Taking  flow  field  results  into  considerations,  the  BC  volume  formations  were  also 
analyzed  numerically  and  shown  in  Figure  64.  The  changing  BC  volume  of  case  26.1 
and  26.2  were  compared  to  base  case  25  and  case  26  (Figure  57).  There  are  two 
distinctly  unexpected  characteristic  observed  for  cases  25  and  26.2  (free)  vs.  26  and  26.1 
(50%  or  wall).  As  mentioned,  case  25  is  base  case  with  free  boundary  and  case  26.1  is 
the  same  except  the  location  of  the  LB.  In  case  25,  the  boundary  is  located  well  outside 
the  BC  max  radius  and  case  26.1  the  boundary  is  located  inside  the  BC  max  radius.  The 
volumetric  formation  of  BC  result  shows  that  results  of  case  25  and  26.2  are  not  identical. 
During  the  initial  collapse  of  BC  near  40  to  50msec,  the  results  show  case  26.2  to 
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collapse  a  lot  faster  than  that  of  ease  25.  Additionally,  near  30msee  and  beyond,  the  ease 
26.2  shows  greater  volume  spike  than  ease  25.  This  is  during  the  seeondary  pressure 
spike  where  BC  hammer  is  present  (approximately  35msec)  and  from  the  observed 
results,  ease  26.1  and  26.2  had  the  most  impaet  in  their  BCs.  For  time  step  ranges  of  50 
to  lOOmsee,  where  LB  reflective  shockwaves  are  present,  ease  26  and  26.1  had  the  most 
impaet.  This  was  expeeted  due  to  presences  of  refleetive  setting  in  their  LB  eonditions. 
Lastly,  time  step  ranges  from  125  to  150mseo  show  oases  26  and  26.1  with  most 
additional  cavitation  formed  due  to  same  reasoning  behind  presence  of  reflective  LB 
eonditions. 


1^  Bulk  Caviation  Volume  for  Case  26  at  Various  Reflectivity 


Figure  64.  Bulk  Cavitation  Volume  for  Various  Lateral  Boundary  Refleotivity 
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VII.  PRECURSOR  TO  FULLY  COUPLED  RUNS 


A.  THE  BLOCKED  CELL  METHOD 

The  presence  of  foreign  body  and  its  effect  on  UNDEX  parameters  in  LOD  is  of 
interest  for  this  chapter.  Likewise,  UNDEX  parameter’s  impact  on  the  ESA  is  also  of 
interest.  As  a  precursor  to  including  rigid  body  or  Lagrangian  solid  in  the  Eulerian 
domain  for  fully  coupled  run,  use  of  “blocked  cell”  method  was  utilized.  The  blocked 
cell  method  is  set  and  invoked  during  Pregemini  phase,  where  inclusion  of  such  object  in 
the  Eulerian  domain  acts  as  a  rigid  body  in  lieu  of  an  actual  Lagrange  structure  in  its 
place,  without  actually  running  a  fully  coupled  run.  The  placements  of  foreign  bodies  in 
three  locations  as  blocked  cells  (BK)  are  shown  in  Figure  65.  Generality  exists  to  show 
that  placements  of  foreign  surface  object  near  origin,  at  midpoint  and  at  the  edge  of  BC 
radius  will  have  an  observable  impact  on  UNDEX  parameters;  hence,  making  these  cases 
worth  the  look  prior  to  fully  coupled  simulations.  As  before,  the  target,  charge  and 
DYSMAS  set  up  are  as  described  in  Table  4. 


A  Y 


Air  Column 


CHARGE 

(X,Y,2)c 


Water  Column 


Figure  65.  Blocked  Cell  Methods 
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Case 

# 

Run 

Mode 

Euler 

Dimesions 
(X,Y,Z)  [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

29 

2D,  Euler 

300, 1,  300 

10 

15,  105,  255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

30 

2D,  Euler 

300, 1,  300 

10 

15 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

31 

2D,  Euler 

300, 1,  300 

10 

105 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

32 

2D,  Euler 

300, 1,  300 

10 

255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Table  4.  Case  Studies  for  Blocked  Cell  Methods 


B,  SUMMARY  OF  RESULTS  AND  ANALYSIS 

1.  Flow  Field  Analysis 

The  flow  field  results  for  cases  29  through  32  are  depicted  in  figure  66  through 
69.  The  base  case,  without  any  blocked  cell  is  shown  in  Figure  66.  No  unusual 
observation  was  made  and  all  time  steps  with  significant  UNDEX  events  following  the 
same  pattern  as  previous  cases  studied:  (a)  initial  SW  propagation  (b)  BC  reaching 
maximum  radius  and  depth  (c)  collapse  of  BC  and  initiation  of  hammer  SWs  (d) 
emergence  of  BB  SW  as  it  travels  to  intercept  water-air  interface  to  invoke  additional 
cavitation.  Two  things  that  stand  out  are  the  characteristics  of  residual  cavitations  shown 
in  Figure  66(d). 


Figure  66.  Case  29  (Base),  Time  Step  (a)  10  (b)  30  (c)  50  (d)  100msec 
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Case  30  with  blocked  cell  closest  to  the  detonation  is  depicted  in  Figure  67.  The 
distinctly  different  characteristics  are  (a)  observation  of  local  cavitation  near  and  around 
the  blocked  cell,  (b)  disruption  in  BC  formation  by  the  blocked  cell,  (c)  presence  of 
additional  hammer  SWs  originating  from  collapse  of  local  cavitation  of  blocked  cell  and 
speed  at  which  initial  BC  collapses. 


Figure  67.  Case  30  (BKl),  Time  Step  (a)  10  (b)  30  (c)  50  (d)  100msec 

Case  31  with  blocked  cell  at  midpoint  to  BC  radius  is  depicted  in  Figure  68.  The 
distinctly  different  characteristics  are  (a)  lack  of  local  cavitation  near  and  around  the 
blocked  cell,  (b)  disruption  in  BC  formation  by  the  blocked  cell  as  it  fully  matures  and 
local  cavitation,  (c)  presence  of  additional  hammer  SWs  originating  from  collapse  of 
local  cavitation  of  blocked  cell  merging  with  the  initial  BC  collapse  and  (d)  emergence  of 
laterally  reflected  SW  originated  from  the  local  cavitation  hammer  shockwave. 
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Figure  68.  Case  31  (BK2),  Time  Step  (a)  10  (b)  30  (c)  50  (d)  100msec 


Case  32  with  blocked  cell  located  well  outside  the  max  radius  of  BC  is  depicted  in 
Figure  69.  Once  again,  the  distinctly  different  characteristics  are  (a)  lack  of  local 
cavitation  near  and  around  the  blocked  cell,  (b)  lack  disruption  in  BC  formation  and 
presence  of  local  cavitation,  (c)  lack  of  additional  hammer  SWs  originating  from  collapse 
of  local  cavitation  and  (d)  emergence  of  laterally  reflected  SW  originated  from  unknown 
source.  The  origin  of  this  laterally  reflected  wave  can  be  speculated  from  local  cavitation 
collapse  from  blocked  cell,  but  with  time  step  and  location  of  blocked  cell  for  the  case  32, 
this  is  unlikely. 


Figure  69.  Case  32  (BK3),  Time  Step  (a)  10  (b)  30  (c)  50  (d)  100msec 
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2,  Target  Pressure  Analysis 

The  target  pressure  profiles  for  eaeh  bloeked  cell  location  are  depicted  in  figures 
70  through  72.  Figure  70  is  target  pressure  profde  for  BKl  compared  to  the  base  case 
without  any  blocked  cells.  As  can  be  seen,  although  there  are  no  impacts  to  the 
magnitude  of  initial  shockwave  of  nearly  700psi  near  10msec,  a  significant  differences  in 
BC  collapse  shockwave’s  time  and  magnitude  can  be  observed.  While  the  base  case’s 
BC  collapsing  shockwave  reaches  close  to  400psi,  the  BKl  only  experiences  close  to 
lOOpsi.  BKl  also  experiences  BC  closure  almost  10msec  faster  than  base  case.  As 
observed  in  the  flow  field  results,  the  presence  of  BKl  also  expedites  the  BC  collapse. 


Pressure  Profile  vs.  Time 


Figure  70.  Pressure  Profile  for  BKl  and  Base  Case 


Target  pressure  profde  for  BK2  and  base  case  is  depicted  in  Figure  71.  Once 
again,  the  initial  shockwaves  are  almost  identical  in  magnitude  and  timing  with  close  to 
240psi  at  time  25msec.  Unlike  previous  case,  this  time  the  BC  collapse  shockwave  is 
higher  for  the  BK2  at  close  to  85psi  while  base  case  hovers  at  25psi.  Post  100msec  show 
presence  of  additional  cavitations  formed  when  shockwaves  originating  from  the  collapse 
of  local  cavitation  is  reflected  off  the  lateral  boundaries,  shown  in  flow  field  of  Figure  68. 
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Pressure  Profile  \s.  Time 


Figure  7 1 .  Pressure  Profile  for  BK2  and  Base  Case 


Pressure  Profile  vs.  Time 


Figure  72.  Pressure  Profile  for  BK3  and  Base  Case 


Target  pressure  profile  for  BK3  and  base  case  is  depicted  in  Figure  72.  Once 
again,  the  initial  shockwaves  are  almost  identical  in  magnitude  and  timing  with  close  to 
lOOpsi  at  time  46msec.  Similar  to  previous  case,  this  time  the  BC  collapse  shockwave  is 

higher  for  the  BK3  at  close  to  35psi  while  base  case  hovers  at  20psi.  Post  120msec  also 

78 


show  presence  of  additional  cavitations  formed  when  shockwaves  originating  from  the 
collapse  of  local  cavitation  is  reflected  off  the  lateral  boundaries,  shown  in  flow  field  of 
Figure  69. 


3,  Vertical  Take-Off  Velocity  Analysis 

The  VTO  for  each  blocked  cell  location  are  depicted  in  figures  73  through  75. 
The  VTO  for  base  case  and  BKl  are  depicted  in  Figure  73.  Both  cases’  VTO  peaks  at 
11  ft/s  followed  by  an  immediate  collapse  and  fluctuation  near  and  around  3ft/s.  The 
VTO  fluctuations  for  the  base  case  seems  to  be  more  dramatic  in  that  lack  of  blocked  cell 
that  can  hinder  movement  will  freely  let  water  to  slush  around  the  domain.  Conversely, 
the  presence  of  blocked  cells  greatly  decreases  the  movement  of  the  water  around  it  since 
the  energy  is  transferred  to  and  from  blocked  cell  invoking  a  beginning  of  FSA 
phenomenon. 


Vertical  Take-Off  Velocity  (VTO)  vs.  Time 


Figure  73.  VTO  for  BKl  and  Base  Case 


The  VTO  for  the  base  case  and  BK2  are  depicted  in  Figure  74.  Both  cases’  VTO 

peaks  at  1.5ft/s  followed  by  an  immediate  collapse  and  fluctuation  near  and  around 

0.25ft/s.  The  VTO  fluctuations  for  the  BK2  seem  to  be  more  dramatic  with  dip  of -1.5ft/s 
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at  approximately  40msec.  The  dip  at  40msec  is  counterintuitive  since  the  blocked  cell  is 
not  moving,  but  further  investigating  flow  field  results  (Figure  68),  show  that  this  is  due 
to  traveling  shockwave  originated  from  the  collapsing  local  cavitation. 


Vertical  Take-Off  Velocity  (VTO)  \s.  Time 


Figure  74.  VTO  for  BK2  and  Base  Case 


Lastly,  VTO  for  base  case  and  BK3  are  compared  in  Figure  75.  For  this  case,  the 
presence  of  initial  shockwaves  at  0.3psi  at  55msec  and  presence  of  shockwave 
originating  from  collapse  of  local  cavitation  are  close  to  being  identical,  except,  for  BK3, 
the  VTO  continues  to  increase  throughout  the  time  step.  Since  blocked  cells  are  treated 
as  a  rigid  material  and  reflecting  surface,  BK3’s  characteristic  can  only  be  explained  as 
error  in  Gemini’s  solution  [12].  Numerous  attempts  were  made  to  fix  this  error  through 
remeshing  Eulerian  domain  and  placement  of  the  blocked  cells,  as  well  as,  target 
position.  At  this  time,  an  acceptable  solution  does  not  exist  although  the  pattern  and 
presence  of  initial  and  additional  shockwave  are  consistent  with  previous  cases. 
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Vertical  Take-Off  Velocity  (VTO)  vs.  Time 


Figure  75.  VTO  for  BK3  and  Base  Case 


4,  Bulk  Cavitation  Analysis 

The  BC  results  for  all  blocked  cell  cases  are  depicted  in  Figure  76.  Similar 
characteristics  to  previous  cases  are  still  observed,  where  BC  is  initiated  at  10msec  when 
initial  shockwave  reaches  water-air  interface,  reaches  its  maximum  radius  at 
approximately  35msec,  dissipates  completely  by  60msec  and  additional  cavitation 
observed  near  120msec  where  BB  reaches  water-air  interface  to  invoke  this  phenomenon. 
Two  distinctly  different  patterns  emerge  for  BK2  where  the  BC  starts  to  collapse  22msec 
with  another  increase  near  40msec  when  local  cavitations  are  invoked  and  completes 
dissipation  near  60msec.  For  BK3,  lack  of  local  cavitation  reduces  the  overall  BC 
volume;  hence,  its  collapse  prior  to  60msec  can  also  be  observed.  All  four  cases 
experienced  additional  cavitation  formations  starting  at  120msec  when  BB  shockwave 
emerged  to  initiate  another  set  of  BC  near  water-air  interface. 
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x  10^  Bulk  Caviation  Volume  vs.  Time 


Figure  76.  Bulk  Cavitation  Volume  for  Blocked  Cell  Cases 
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VIII.  CONCLUSION  AND  RECOMMENDATION  FOR  FURTHER 

STUDIES 


A,  CONCLUSIONS 

An  in-depth  analysis  and  characterization  of  littoral  ocean  domain  (LOD)  were 
conducted  for  transient  phase  of  UNDEX  by  analyzing  resulting  parameters  such  as 
initial  shock  wave,  vertical  take-off  velocity  and  bulk  cavitation.  Implementing  various 
conditions  such  as  ocean  depth,  charge  size  and  depth,  boundaries  for  bottom  and  side 
conditions  of  fluid  domain  and  even  including  a  small  rigid  body,  these  results  were 
compared  and  contrasted  to  the  current  Eulerian  fluid  model. 

Eirst,  results  of  varying  ocean  depths  show  clear  distinction  in  UNDEX 
characterization  at  depth  shallower  than  300ft.  At  this  depth,  multiple  shockwaves  and 
resulting  vertical  take-off  velocity  exists  due  to  additionally  induced  cavitation  resulting 
from  the  reflecting  pressure  wave  from  the  shallower  bottom.  The  reflecting  shockwave 
and  energy  contained  within  are  the  main  culprit  of  shallow  UNDEX  characteristics. 

Second,  the  results  of  varying  charge  size  and  depth  showed  that  charge  size  of 
less  than  or  equal  to  3001bs  showed  linear  relationships  and  as  the  charge  depth  reached 
water-air  or  water-bottom  interface,  the  characteristics  of  resulting  UNDEX  parameters 
became  increasingly  amplified  and  chaotic.  The  combined  bottom  bounce  and  initial 
shockwave  for  charge  depth  at  water-bottom  interface  showed  combined  effects  of  its 
resulting  UNDEX  parameters  with  greater  BC  volume  while  delayed  initial  response. 

Third,  results  of  varying  lateral  boundary  conditions  show  that  as  the  lateral 
boundary  distance  is  brought  closer  to  the  detonation  source,  inside  the  radius  of  bulk 
cavitation,  its  behavior  and  the  rest  of  the  UNDEX  parameters  also  become  increasingly 
amplified  due  to  similar  effects  seen  in  shallower  bottom  depths.  Varying  its  reflectivity 
also  amplified  UNDEX  reactions;  however,  differences  in  reflectivity  between  free  and 
50%  lateral  boundary  condition  inside  and  outside  the  radius  of  BC  showed  negligible 
changes  in  the  UNDEX  results. 
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Lastly,  adding  blocked  cells  prior  to  a  full  scale  coupled  run  showed  that  fluid 
behaves  more  erratically  as  these  blocked  cells  are  situated  within  the  radius  of  forming 
BC.  As  fluid  behaves  differently  around  the  existing  small  rigid  blocks,  so  too  are  the 
UNDEX  parameters  due  to  reflective  energy  from  these  rigid  blocks,  invoking  energy 
back  into  the  fluid  domain.  As  a  result,  a  precursor  to  how  fluid  may  behave  starts  to 
emerge  for  a  fluid  structure  interaction  or  fully  coupled  run  cases. 

As  a  result,  recommended  characterization  and  boundaries  of  LOD  for  UNDEX 
studies  start  to  emerge.  Erom  this  study,  the  ocean  depth  of  300ft,  lateral  boundary 
distance  of  greater  than  1/3  radius  distance  from  the  edge  of  bulk  cavitation,  charge 
weight  of  3001bs  or  less  and  its  placement  at  minimum  of  1/3  distance  away  from  the 
existing  boundaries,  such  as  bottom  or  lateral  wall  are  recommended.  Eor  fully  coupled 
runs,  the  placement  of  Eagrangian  solid  within  the  bulk  cavitation  radius  and  inclusion  of 
various  bottom  sediments  to  increase  realism  with  caveats  of  setting  its  reflectivity 
condition  to  rigid  or  wall  are  also  recommend  for  EOD.  Although  cyclic  bubble  was  not 
studied  extensively  during  this  research,  the  obvious  choice  of  bubble  setting  is  to  pick  a 
charge  size  that  will  produce  an  initial  bubble  size  consistent  with  previous  EOD 
boundaries.  If  purpose  is  to  study  the  transient  phase  of  UNDEX  then  this  is  not  a  huge 
factor.  However,  if  whipping  analysis  is  of  interest,  choosing  a  charge  size  that  will 
invoke  a  bubble  radius  of  at  minimum  of  less  than  chosen  EOD  depth  is  recommended. 

B,  RECOMMENDATION  FOR  FURTHER  STUDIES 

While  this  study  shed  some  light  on  the  benefits  of  establishing  the  boundaries  of 
EOD  for  future  coupled  runs  of  various  Eagrangian  models  or  an  actual  ECS  EE  model, 
due  to  the  sheer  size  of  the  problem  in  computing  time  and  efficiency,  other  parts  of 
Eulerian  characteristics  such  as  implementing  viscous  code,  adding  shelves  within  the 
fluid  domain  to  study  the  effects  of  its  angle  and  further  study  of  steady  state  UNDEX 
phenomena  like  cyclic  bubbles  and  its  effect  are  recommended. 
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APPENDIX  A.  COMBINED  DYSMAS/GEMINI  SIMULATED  CASES 


Case 

# 

Run 

Mode 

Euler 

Dimesions 

(X,Y,Z)  [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

Remarks 

1 

2D,  Euler 

800,  1,  1000 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Varying  Ocean  Depth  from  1000ft  to  75ft 
deep.  All  other  UNDEX  parameters, 
including  charge  size  and  charge  depth 

maintained  constant. 

2 

2D,  Euler 

800,  1,  900 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

3 

2D,  Euler 

800,  1,  800 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

4 

2D,  Euler 

800,  1,  700 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

5 

2D,  Euler 

800,  1,  600 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

6 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

7 

2D,  Euler 

800,  1,  400 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

8 

2D,  Euler 

800,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

9 

2D,  Euler 

800,  1,  200 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

10 

2D,  Euler 

800,  1,  175 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

11 

2D,  Euler 

800,  1,  150 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

12 

2D,  Euler 

800,  1,  125 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

13 

2D,  Euler 

800,  1,  100 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

14 

2D,  Euler 

800,  1,  75 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

15 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Varying  Charge  Size.  All  other  UNDEX 
parameters  maintained  constant. 

16 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

200 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

17 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

300 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

18 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

400 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

19 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

500 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

20 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

100 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Varying  Charge  Depth.  All  other  UNDEX 
parameters  maintained  constant. 

21 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

200 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

22 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

300 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

23 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

400 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

24 

2D,  Euler 

800,  1,  500 

0.820 

100 

HBX-1 

100 

500 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

25 

2D,  Euler 

230,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Base  Case,  free  boundary  @  xmax 

26 

2D,  Euler 

164,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

Varying  Case,  wall  boundary  @  xmax.  For 
26.1  &  26.2,  varied  reflectivity  50%  &free 
for  boundary  located  inside  BC  zone. 

26.1 

2D,  Euler 

164,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

0.50 

Wall 

Wall 

0.656 

26.2 

2D,  Euler 

164,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

27 

2D,  Euler 

197,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

28 

2D,  Euler 

230,  1,  300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

29 

2D,  Euler 

300,  1,  300 

10 

15,  105,  255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Base  Case,  No  Blocked  Cell 

30 

2D,  Euler 

300,  1,  300 

10 

15 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Blocked  cell  at  position  (1) 

31 

2D,  Euler 

300,  1,  300 

10 

105 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Blocked  cell  at  position  (2) 

32 

2D,  Euler 

300,  1,  300 

10 

255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Blocked  cell  at  position  (3) 
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APPENDIX  B:  BULK  CAVITATION  THEORY 


The  following  section  describes  how  the  method  of  Arons  is  used  to  map  out  the 
bounds  of  the  cavitated  region  that  forms  when  surface  cutoff  tries  to  lower  the  absolute 
pressure  below  the  cavitation  pressure.  A  continuation  from  section  II.  B.  2  of  this  study, 
it  is  a  direct  summary  from  the  critical  sections  within  the  technical  paper  published  by 
Costanzo  and  Gordon  in  May  1983,  titled  “A  Solution  to  the  Axisymmetric  Bulk 
Cavitation  Problem”  and  for  that  I  am  grateful  for  its  use  and  many  thanks  to  Contanzo 
and  Gordon  for  use  of  their  work  as  primary  reference  in  characterizing  BC  for  this  study 

[9]. 

A,  DERIVATION  OF  METHODS  OF  ARONS  AND  TANGENT  RULE 

The  target  pressure  (P)  at  its  location  of  (X,  Y)  is  assumed  to  be  the  absolute 
pressure  prior  to  arrival  of  the  reflected  wave  from  water-air  interface.  This  pressure  is 
culmination  of  overpressure  generated  by  SW,  atmospheric  and  hydrostatic  pressure. 
Now,  let  P  be  function  of  radial  (r)  and  angular  (a)  coordinate  originating  from  the  image 
charge,  Wi,  as  shown  in  Figure  20.  The  upper  cavitation  boundary  is  defined  as  the  locus 
of  points  at  which  the  absolute  pressure  drops  to  the  cavitation  pressure  upon  arrival  of 
the  rarefaction  wave  [9].  At  this  point,  water  at  cavitation  will  remain  cavitated  as  long 
as  the  absolute  pressure  does  not  rise  above  the  vapor  pressure  of  water.  For  all  practical 
purpose  for  the  method  of  Arons,  the  vapor  and  cavitation  pressure  are  assumed  to  be 
zero.  Thus,  the  equation  defining  the  upper  cavitation  boundary  is 


P[a,r)  +  P^  =0 


(B.l) 


If  Pr,  the  reflected  wave,  is  expressed  in  terms  of  the  charge  weight  and  standoff 

1/3  B 

using  the  method  of  images,  then  Pr  =  -A(W  /r)  .  Substituting  this  expression  into 
Equation  (B.l)  gives 


P[a,r^-A 


=  0 


V  r  j 
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(B.2) 


Where,  A  and  B  are  constants  specific  to  the  charge  or  explosive  types.  Now,  before  the 
cavitation  boundary  is  discussed,  the  term  breaking  pressure  must  be  defined.  The 
breaking  pressure  is  the  magnitude  of  the  rarefaction  wave  (or  relief  wave,  as  commonly 
termed)  which  reduces  the  absolute  pressure  at  a  point  to  the  cavitation  pressure.  In  other 
words,  since  the  cavitation  pressure  is  taken  to  be  0  psi  absolute,  the  breaking  pressure 
has  a  magnitude  equal  to  the  absolute  pressure  at  a  point  prior  to  the  occurrence  of 
cavitation  at  that  point.  The  equation  of  the  lower  cavitation  boundary  is  derived  from 
consideration  of  the  propagation  of  this  breaking  pressure  into  the  uncavitated  water 
beneath  this  boundary.  Let  Pr  =  P(a,  r)  be  the  breaking  pressure  for  a  point  lying  at  the 
lower  cavitation  boundary.  At  this  lower  boundary,  P(a,  r)  must  propagate  into 
uncavitated  water  with  spherical  attenuation  resulting  in  Pr  =  P(a,  r)(R/r)®.  This  is 
represented  in  Figure  77. 


Figure  77.  Propagation  of  Breaking  Pressure  into  Uncavitated  Water,  from  [9] 

This  pressure  expression  must  have  a  faster  decay  rate  than  the  general  expression 
for  absolute  pressure,  P  =  P(a,  r),  or  the  water  located  along  the  dashed  line  extending 
below  the  lower  cavitation  boundary  will  continue  to  cavitate.  Thus,  at  a  point  on  the 
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lower  cavitation  boundary,  the  decay  rate  of  both  of  these  pressure  expressions  is  the 
same.  Therefore: 


dr' 


(B.2.a) 


If  all  terms  of  this  equation  are  shifted  to  one  side  of  the  equal  sign  and  r‘®  is 
factored  out,  then  application  of  the  chain  rule  to  differentiation  of  a  product  results  in: 


d 

dr 


[r^)p{a,r)-P{a,r){R^) 


\r=R 


0 


(B.2.b) 


This  expression  may  be  simplified  further  by  recognizing  that  the  product, 
P(a,r)(R®),  pertains  to  a  specific  point  and  thus  may  be  treated  as  a  constant  when 
differentiating.  Therefore,  this  equation  becomes: 


=  0 

(B.3) 

where  r  =  R,  the  point  at  the  lower  cavitation  boundary  for  a  given  a.  Equation  (B.3)  is 
the  method  of  Arons  for  determining  the  lower  cavitation  boundary. 

In  Figure  78,  the  general  shape  of  the  BC  envelope  as  determined  by  the  method 
of  Aron  is  given.  It  is  of  fundamental  interest  to  the  BC  problem  to  locate  the  point  at 
which  a  line  drawn  from  the  image  charge  is  tangent  to  the  upper  cavitation  boundary. 
For  the  derivation  of  this  tangent  point,  both  sides  of  Equation  (B.2)  are  multiplied  by  r® 
and  the  resulting  equation  and  its  total  derivation  with  respect  to  r  are  written  below  as 
equations  (B.4)  and  (B.5),  respectively. 


—  (r^)P{a,r) 


[r^)p{a,r)-A[w^'^)'' 


(B.4) 


A. 

dr  L 


[r'')p{a,r)- 


■  _  5 

\ir^\pia,r) 

d 

+ — 

\ir^\pia,r) 

^  da^ 

dr 

LV  /  V  ’  /J 

da 

LV  /  V  /J 

V  Sr  y 

=  0 


(B.5) 
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For  a  given  value  of  da/dr,  the  simultaneous  solution  of  these  two  equations  gives 
the  eorresponding  values  of  a  and  r.  When  the  ray  extending  from  the  image  eharge  is 
tangent  to  the  upper  eavitation  boundary,  a  is  a  maximum  and  da/dr  =  0.  When  this 
value  for  da/dr  is  substituted  into  Equation  (B.5),  the  following  is  obtained: 


dr 


(r’^)P{a,r) 


0 


(B.6) 


Therefore,  the  simultaneous  solutions  of  Equations  (B.4)  and  (B.6)  give  the  value 
of  a  and  r  at  the  tangent  point.  Equation  (B.6)  also  happens  to  be  equivalent  to  Equation 
(B.3),  the  equation  whose  solution  determines  the  lower  eavitation  boundary.  This 
indieates  that  the  lower  and  upper  eavitation  boundaries  interseet  at  the  point  at  whieh  a 
line  drawn  from  the  image  eharge  is  tangent  to  the  upper  boundary.  This  rule  of  tangeney 
is  also  illustrated  in  Eigure  78. 


Eigure  78.  Bulk  Cavitation  Bounds  and  Rule  of  Tangeney,  from  [10] 

B.  DEVELOPMENT  OF  THE  CLOSURE  MODEL 
1.  General  Description 

The  general  representation  of  a  point  whieh  lies  within  the  eavitated  region  is 
given  in  Figure  79.  Upon  the  arrival  of  the  relief  wave  at  this  point,  the  vertieal 


eomponent  of  the  instantaneous  water  partiele  veloeity  is  the  veetor  sum  of  the  vertieal 
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components  of  the  two  veloeity  veetors  as  shown.  This  veetor  sum  is  termed  the  vertieal 
water  partiele  kiekoff  veloeity  or  vertieal  take-off  veloeity  (VTO)  at  point  (X,  Y)  and  is 
dependent  solely  upon  the  magnitude  and  directions  of  the  ineident  and  refleeted  acoustic 
water  partiele  veloeities.  The  direetions  of  these  veloeities  are  defined  by  the  unit 

veetors,  ^  and  J ,  as  shown  in  the  figure. 


At  every  horizontal  range  within  the  extent  of  the  BC  envelope,  there  exists  a 
column  of  water  whieh  lies  between  the  surface  and  the  upper  cavitation  boundary.  Since 
this  water  does  not  eavitate,  it  is  assumed  that  all  water  partieles  eontained  in  this  vertical 
column  are  kieked  off  simultaneously  with  kiekoff  velocity  of  the  water  partiele  located 
at  the  upper  eavitation  boundary.  This  is  illustrated  in  Figure  80.  This  surfaee  layer, 
whieh  is  kieked  off  at  the  time  of  relief  wave  arrival  at  the  upper  eavitation  boundary 

with  an  initial  veloeity,  ,  is  deeelerated  by  atmospherie  pressure  and  gravity.  At  a  short 
time  later,  the  water  partiele  lying  direetly  beneath  this  surfaee  layer  is  kieked  off  with  an 

initial  veloeity,  . 

Sinee  this  partiele  lies  within  the  eavitated  region  and  is  kieked  off  after  the 

surfaee  layer,  it  beeomes  separated  from  the  surfaee  layer  and  therefore  has  no 

atmospherie  pressure  aeting  upon  it.  Thus,  only  gravity  deeelerates  this  partiele. 

Eventually,  this  water  partiele  will  eollide  with  the  surface  layer  above  it.  In  the 
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development  of  the  closure  model,  it  is  assumed  that  this  is  a  perfectly  inelastic  collision. 
Thus,  this  particle  and  the  surface  layer  above  it  now  form  an  augmented  surface  layer 
which  has  a  velocity  derived  from  the  conservation  of  momentum  consideration.  As  with 
the  original  surface  layer,  this  augmented  surface  layer  has  atmospheric  pressure  and 
gravity  acting  to  decelerate  it. 


Figure  80.  Kickoff  Velocity  of  the  Surface  Layer,  from  [9] 

Since  the  particles  lying  below  the  original  surface  layer  are  all  kicked  off  at 
different  times  with  different  vertical  kickoff  velocities,  inelastic  collisions  will  occur  one 
at  a  time  between  the  growing  surface  layer  and  the  particles  directly  below  it.  This 
process  is  known  as  accretion.  If  the  surface  layer  displacement  history  at  a  horizontal 
range,  X,  is  plotted  with  t=0  referring  to  the  time  of  explosives  charge  detonation,  the 
curve  in  Figure  81  is  obtained.  This  curve  is  not  quite  a  perfect  parabola  due  to  the  fact 
that  the  surface  layer  mass  is  changing.  Also  note  that  this  curve  accounts  for  the  initial 
displacement  of  the  surface  layer  due  to  the  incident  SW  prior  to  relief  wave  arrival  at  the 
upper  cavitation  boundary. 
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Figure  81.  Surface  layer  Displacement,  from  [9] 


The  vertical  component  of  the  water  particle  velocity  for  the  point  which  lies  at 
the  lower  cavitation  boundary  at  this  same  horizontal  range,  X,  also  must  be  determined. 
Since  cavitation  does  not  extend  below  this  point,  separation  will  not  occur  between  the 
underlying  water  particles.  Hence,  the  water  which  lies  below  this  point  at  the  lower 
cavitation  boundary  will  have  a  vertical  velocity  component  which  is  dependent  upon  the 

varying  velocities  and  ^  and  their  corresponding  afterflow  terms.  This  is  depicted  for 
point  (X,  Yl)  in  Figure  82. 


Figure  82.  Velocity  of  a  Point  at  the  Lower  Cavitation  Boundary,  from  [9] 
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2,  Bulk  Cavitation  Closure  Pulse 

The  magnitude  of  the  pressure  pulse  generated  by  the  collision  of  the  accreting 
lower  cavitation  boundary  at  closure  will  now  be  determined.  Since  in  the  closure  model, 
the  vertical  components  of  the  velocities  of  both  the  surface  layer  and  the  lower 
cavitation  boundary  for  a  particular  horizontal  range,  X,  are  known  at  the  time  of  closure, 
the  closure  pressure  pulse  magnitude.  Pc,  can  be  calculated  by  multiplying  one  half  the 
relative  velocity  of  impact  of  these  two  water  masses  by  the  characteristic  impedance  of 


the  medium.  That  is: 


(B.7) 


where  V  is  positive  upward  and  Vu  and  Vl  refer  to  the  vertical  velocity  components  at 
closure  of  the  surface  or  upper  layer  and  the  lower  cavitation  boundary,  respectively.  The 
closure  pressure  pulse  consists  of  two  compressive  waves  of  magnitude.  Pc;  one  that 
travels  upward  and  one  that  travels  downwards.  These  closure  pressures  are  most  readily 
superimposed  upon  the  incident  SW  and  reflected  pressure  at  the  horizontal  range  at 
which  closure  first  occurs,  since  at  the  range,  the  direction  of  propagation  of  the  closure 
pulse  is  vertical.  At  other  horizontal  ranges,  the  direction  of  propagation  of  the  closure 
pulse  varies  from  the  vertical  direction  by  an  angle  \|/,  as  shown  in  Figure  83.  According 
to  the  work  done  by  Cushing  [10],  if  C  is  the  sound  speed  in  water  and  Vc  is  the  speed  at 
which  the  cavitated  region  is  closing  at  some  horizontal  range,  X,  then: 


smi//  = 


(B.8) 


Thus,  the  magnitude  of  the  pressure  pulse  which  propagates  along  a  path  oriented 
at  this  angle  \\i  is  determined  by  dividing  the  expression  of  Equation  (B.7)  by  cos  \\i.  That 


pC(V,-Vu) 

2  cosy/ 


(B.9) 
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Figure  83.  Closure  Pulse  at  a  Horizontal  Range  Different  from  that  of  First 

Closure,  from  [9] 


The  duration  time  of  the  elosure  pulse  at  the  depth  of  elosure  is  determined  by  the 
method  of  images  and  is  expressed  by: 


-cost//' 


(B.IO) 


where  Dc  is  the  depth  of  elosure  at  the  horizontal  range  of  interest,  X.  Thus,  regardless 
of  the  angle  of  propagation  of  the  closure  pulse,  the  impulse  (pressure  X  duration) 
associated  with  cavitation  closure  at  a  particular  horizontal  range  is  the  same  and  may  be 
expressed  as  shown  in  Equation  (B.IO). 


Im pulse  -  p{y^  -  Vjj)D^ 


(B.ll) 
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APPENDIX  C.  BULK  CAVITATION  CASE  STUDIES 


Following  case  studies  were  conducted  in  order  to  better  characterize  behaviors  of 
bulk  cavitation  prior  to  its  inclusion  into  shallow  water  or  littoral  ocean  domain 
environments.  The  MATLAB  codes  used  for  this  study  is  included  in  the  Appendix  of 
this  report. 

1,  Comparison  of  2-D  BC  Zone,  varying  charge  type,  charge  weight  and  charge 
detonation  depth  (Appendix  C  contains  MATLAB  codes  for  2D  BC  Zones): 

a,  200  Ih  charges  of  HBX-1,  TNT  and  PETN  at  a  depth  of  25  ft: 

Figure  84  is  a  2D  representation  of  the  bulk  cavitaion  zone  of  a  2001bs  charge 
denotated  at  a  depth  of  25ft.  The  charge  lies  along  the  left  vertical  edge  of  each  subplot. 
Each  subplot  represents  a  different  explosive  type  as  titled.  This  figure  was  created  for 
the  sole  purpose  of  comparing  the  BC  zone  of  the  three  different  types  of  explosive  with 
the  other  variables  (depth  of  charge,  weight  of  charge)  held  constant. 


Figure  84.  2-D  Bulk  Cavitation  of  2001bs  of  FlBX-1,  TNT  and  PETN  @  25ft 

Since  water  cannot  support  a  tension  (i.e.  negative  pressure),  cavitation  occurs  at 
points  in  the  water  column  where  the  sum  of  the  pressures  at  that  point  are  less  than  or 
equal  to  zero.  This  pressure  is  sum  of  atmospheric,  water  depth,  incident  shock  pressure, 
and  rarefaction  wave  pressure.  As  can  be  clearly  seen  from  the  above  figure,  even 


97 


though  all  the  explosive  charge  weights  and  detonation  depths  are  the  same,  each 
explosive  charge  type  creates  a  different  BC  region  in  both  depth  and  radial  distance 
outward.  This  is  due  to  the  fact  that  they  have  different  constants  (Ki,  K2,  Ai,  A2) 
associated  with  them  that  are  used  in  the  calculation  of  their  maximum  pressure  and  their 
decay  constant. 

The  Pentolite  (PETN)  ultimately  had  the  lowest  K2  and  A2  thus  causing  it  to  have 
the  deepest  but  shortest  BC  zone.  TNT  and  HBX-1  are  very  similar  in  their  K2  value 
however  TNT’s  A2  value  is  significantly  smaller  than  HBX-1,  causing  it  to  “reach” 
further  outward  from  point  of  detonation.  In  addition,  the  difference  in  Ai  and  A2  values 
seems  to  be  the  major  driving  factor  for  how  far  out  radially  the  BC  zone  expands.  The 
Pentolite  had  the  largest  difference  and  thus  the  smallest  radial  distance,  while  TNT  has 
the  smallest  difference  and  thus  the  largest  radial  distance  of  the  BC  zone. 

b.  100  lb,  200  lb,  and  300  lb  charges  of  HBX-1  at  50  ft: 

Figure  85  is  a  2D  representation  of  the  bulk  cavitaion  zone  of  a  HBX  explosive 
charge  of  varying  weights  (1001b,  200  lb,  and  3001b)  detonated  at  a  constant  depth  of  50 
feet.  Once  again  the  location  of  explosive  depicted  in  the  left  vertical  axis  of  each 
subplot.  The  purpose  of  the  plot  is  to  represent  the  effect  of  charge  weight  on  the  bulk 
caviation  zone  with  other  variables  (charge  type,  detonation  depth)  held  constant.  As 
expected,  as  the  charge  weight  is  increased  while  the  detonation  depth  and  charge  type 
are  held  constant,  both  the  depth  and  breadth  of  the  BC  zone  are  increased. 


Figure  85.  2-D  Bulk  Cavitation  of  100,  200,  and  3001bs  of  HBX-1@  50ft 
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c,  250  lb  charge  of  TNT  detonated  at  varying  depths  (  5,  50  and  500ft): 


Figure  86  is  a  2D  representation  of  the  bulk  cavitaion  zone  of  a  TNT  explosive 
charge  of  2501bs  detonated  at  varying  depths  of  5,  50,  and  500  ft.  Figure  87  shows 
extended  range  of  250  lb  at  500ft  for  obvious  reason.  Once  again,  the  explosive  is 
depicted  along  the  left  vertical  axis  of  each  subplot.  Again,  the  purpose  of  the  plot  is  to 
represent  the  effect  of  charge  weight  on  the  bulk  caviation  zone  with  other  variables 
(charge  type,  detonation  depth)  held  constant. 
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Figure  86.  2-D  Bulk  Cavitation  of  2501bs  TNT  detonated  @  5,  50  and  500ft 
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Figure  87.  Extended  View  of  2501bs  TNT  detonated  @  500ft 


As  expected,  as  the  depth  of  the  detonation  point  of  a  charge  varies,  with  the 
weight  and  type  of  charge  being  held  constant,  the  BC  region  grows  both  in  depth  and 
breadth,  with  the  exception  of  the  500’  depth  charge.  This  charge  does  have  massive 
radial  distance  (more  than  twice  the  radial  distance  of  charge  at  50ft)  compared  to  the 
other  two  and  the  thickness  of  the  BC  zone  seems  to  grow  radially  rather  than  maximum 
thickness  being  at  or  near  the  center  of  the  BC  zone.  This  is  most  likely  due  to  the  fact 
that  the  pressure  from  incident  and  image  charges  are  very  low  (having  decayed 
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significantly  already  as  it  traveled  vertieally  500ft)  and  are  at  a  relatively  shallow  depth. 
Therefore,  the  total  pressure  at  this  water  eolumn,  although  significantly  lower  than  other 
two  depth  cases,  is  enough  to  ereate  a  eavitation  while  the  surrounding  pressure  rises  to 
the  vapor  pressure  of  the  water  (about  0.3  psi).  In  addition,  the  SW  that  causes  the 
pressure  havoc  in  this  water  column  and  causes  the  BC  has  traveled  significantly  farther 
than  shallower  depth  eharges. 

2,  3-D  visualization  of  the  BC  space  for  an  underwater  explosion  resulting  from  the 
detonation  of  a  l,0001hs  charge  of  TNT  located  at  a  depth  of  25ft:  (Appendix  C 
contains  MATLAB  codes  for  3D  TNT  BC  Zone): 

Figures  88  and  89  depict  the  2D  and  3D  representation  of  l,0001bs  charge  of  TNT 
located  at  25ft  depth.  The  growth  and  shape  of  the  BC  zone  follows  the  same  pattern 
observed  previously  in  section  1  of  this  report.  Figure  5  below  shows  regional  and  line 
depietion  of  the  BC  region.  Following  the  same  trend  of  BC  formation,  the  maximum 
thiekness  appears  at  200ft  radial  distance  from  “ground  zero”  or  point  of  detonation  with 
about  40ft  thickness.  The  radial  length  of  the  BC  region  is  about  650ft. 
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Figure  88.  SPY  and  Line  View  of  l,0001bs  TNT  at  25ft  Depth 

Figure  89  is  a  3D  representation  of  BC.  The  overall  shape  is  almost  like  a  flying 
saucer  with  a  dip  on  the  bottom  side  of  the  dise.  For  the  below  representations,  the  top  of 
the  BC  region  was  removed  to  show  the  eontour  shape  within  the  region.  The  figure  on 
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the  right  depicts  “ground  zero”  and  shows  the  radial  expansion  of  the  explosion  below  the 
BC  zone.  The  usual  maximum  dip  that  is  represented  at  the  center  of  the  BC  zone  is  also 
depicted  in  both  figures  almost  like  an  upside  down  bell. 


Figure  89.  3D  General  Top  and  Peak  Top  View  of  lOOOlbs  TNT  at  25ft  Depth 
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APPENDIX  D.  BULK  CAVITATION  MATLAB  CODES 


Computation  of  Bulk  Cavitation  Zone 
for  Underwater  Explosions 
By 

LT  Dean  Ahn,  USN 


%  This  program  is  used  to  develop  the  BC  envelope  for  an 

%  underwater  charge  of  TNT  at  various  1000  Ibf  charge  weights,  25ft  depth 
%  and  plot  them  in  2D  &  3D. 


%%  Prep  Matlab  workspace 
clear  all; 
clc;  close  all ; 

%  2-D  SPY:  BC  visualization  of  UNDEX,  1000  lb  charge  TNT  @25ft 

clear  all; 

clc; 

%  Constants 

Pa  =  14.7;  %  Atmpospheric  pressure  (psi) 

Gamma  =  62.5/144;  %  Weight  density  of  water  (lb/ft^3) 

C  =  5;  %  Acoustic  velocity  of  water  (ft/msec) 

K1  =  22505;  %  Pmax 
A1  =  1.18;  %  Pmax 
K2  =  0.058;  %  Decay  constant 
A2  =  -0.185;  %  Decay  constant 

%  Plotting  BC  in  2D,  TNT,  1000  lb  @  25ft; 
counter  =  0; 
i  =  1; 

for  W  =  [1000];  %  Equivalent  charge  weights  (lb) 
for  Dl=  [25];  %  Charge  location  depths  (ft) 
counter  =  counter+1; 

A=zeros (50,700) ; 
for  y  =  1:51; 

for  X  =  1:701; 

R1  =  sqrt((Dl  -  (y-l))^2  +  (x-l)^2); 

%  Distance  from  charge  to  desired  location  (ft) 

R2  =  sqrt((Dl  +  (y-l))''2  +  (x-l)^2); 

%  Distance  from  image  charge  to  desired  location  (ft) 
theta  =  K2* (W^ (1/3) ) * ( ( (W^ (1/3) ) /Rl) ^  (A2) )  ; 

%  Decay  Constant  (msec) 

Pi  =  (Kl* (W^ (1/3) /Rl) ^ (Al) ) * (exp  (- (R2  -R1 ) / (C*theta ) ) ) ; 

%  Incident  Pressure  Wave  (psi) 

Ph  =  Gamma* (y-1);  %  Hydrostatic  Pressure  at  y  (psi) 

Pr  =  (Kl*  (  (W^  ( 1 /3 ) /R2 ) '^  (A1 )  )  )  ;  %  Refelcted  Pressure  Wave  (psi) 

F  =  Pi  +  Pa  +  Ph  -  Pr;  %  Upper  Bulk  Caviataion  Boundary 

G1  =  -Pi/  (C*theta) *(!+((  (R2-2*D1* (  (D1+  (y-1) ) /R2)  ) /Rl) * (A2*R2/R1-A2- 

1 ) )  )  ; 

G2  =  - (Al*Pi/R1^2) * (R2-2*D1* (  (D1+ (y-1) ) /R2)  ) ; 

G3  =  Gamma* ( (D1+ (y-1) ) /R2)  ; 

G4  =  (A1/R2) * (Pi+Pa  +  Ph); 

G  =  G1  +  G2  +  G3  +  G4 ;  %  Lower  BC  Boundary 
if  F  >  0.001;  %  Combine  BC  Boundaries 
if  G  <  0; 

_ A(y,x)  =  1; _ 
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end 

end 

if  G  >  0; 

A(y,x)  =  1; 

end 

end 

end 

temp counter )  =  A; 
end 

charge=num2 str (W) ; 
i=i+l; 

end 

f igure4=f igure ( ' NumberTitle ' ,  'off',  'Name',  'Bulk  Cavitation 

lOOOLB  TNT  Charge  UNDEX  at  25FT  Depth',  ' PaperOrientation ' ,  'Landscape') 

subplot l=subplot (2,1,1, ' Parent ' , figure4 ) ; 

box  (subplotl,  'on' ) ; 

hold ( subplotl , ' all ' ) ; 

spy (temp (;,:,!)) 

title ('SPY  View') 

xlabel (' Radius  (ft)') 

ylabel (' Depth  (ft)') 

plot ( figure4 ) 

%%  2-D  Line 

gamma=0 . 037 03 ;  %  seawater  weight  density  (lb/in^3) 

pa=14.7;  %  atmospheric  pressure  (psi) 

c=5;  %  acoustic  velocity  (ft/msec) 

charge_type  =  2 ; 

D  =  2  5; 

W  =  1000; 

%  TNT  Data: 

chg_name= ' TNT  charge ' ; 
depth=num2str (D) ; 
weight=num2 str (W) ; 

Kl=22505; 

Al=1.18; 

K2=0. 058; 

A2=-0 .185; 
ub_data= [ ] ; 
lb_data= [ ] ; 

%  Calculation  of  the  upper  boundary; 
for  x=0 : 700 

for  y=0:0.1:510 

rl=sqrt(  (D-y)''2+x^2)  ; 
r2=sqrt ( (D+y) ^2+x^2 ) ; 
theta=K2*W^ (1/3) * (W^ (1/3) /rl) ^A2  ; 

F1=K1* (W^ (1/3) /rl) ^Al*exp (- (r2-rl) / (c*theta) ) ; 

F2= (gamma*y*12 ) -Kl*(W^(l/3)/r2)^Al+pa  ; 

F=F1+F2  ; 
if  F<=0 

ub_data= [ub_data; F  x  -y]  ; 
break 

end 

end 

end 

%  Calculation  of  the  lower  boundary; 
for  x=0 : 700 


Region : 
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for  y=0:0.1:510 

rl  =  sqrt  (  (D-y)  ''2+x''2  )  ; 
r2  =  sqrt  (  (D+y)  ''2+x''2  )  ; 
theta=K2*W^ (1/3) * (W^ (1/3) /rl) ^A2; 
pi=Kl* (W^ (1/3) /rl) ^Al*exp (- (r2-rl) / (c* theta) ) ; 

Gl=- (pi/ (c*theta) )*(!+(( (r2- (2*D* (D+y) /r2 ) ) /rl ) * ( ( (A2*r2 ) /rl ) -A2-1 ) ) ) ; 
G2=- ( (Al*pi) /rl^2) * (r2-2*D* ( (D+y) /r2) ) ; 

G3= (gamma* 12 ) * ( ( D+y ) /r2 ) ; 

G4=(Al/r2) * (pi+pa+ (gamma*y*12 ) ) ; 

G=G1+G2+G3+G4; 
if  G>=0 

lb_data= [lb_data;G  x  -y]  ; 
break 

end 

end 

end 

%  Plot  upper  and  lower  boundary  data  in  2D: 
subplot  (2,1,2) 

plot (ub_data ( ; , 2 ) , ub_data ( ; , 3 ) , 'b ' , lb_data ( : , 2 ) , lb_data ( : , 3 ) , ' r ' ) 

xlabel (' Radius  [ft]') 

ylabel (' Depth  [ft]') 

title (['2D  Line  View']) 

grid  on 

%%  3D  Mesh  &  Peak  Mesh:  BC  View  of  1000  Ibf  TNT  Charge  at  25ft 
figure 

r  =  ub_data  ( :  , 2  )  ; 

z  =  lb_data  ( : , 3 )  ; 

w  =  ub_data  ( : , 3 )  ; 

theta  =  0 :pi/20 : 2*pi; 

X  =  bsxfun (@times, r, cos (theta) ) ; 

Y  =  bsxfun (Etimes, r, sin  (theta) ) ; 

Z  =  repmat (z, 1, length (theta) ) ; 

W  =  repmat (w, 1, length (theta) ) ; 
grid  on 

%  subplot (1,2,1) 
mesh(X,  Y,  Z,  W) 
title  (  ' 3D  General  View  ' ) 
xlabel ( 'Radius  (ft)') 
ylabel (' Radius  (ft)') 
zlabel  ( 'Depth  [ft]') 

%  subplot (1,2,2) 
figure 

meshc (X,  Y,  Z,  W) 
title (' 3D  Peak  View') 

xlabel ( 'Radius  (ft)') 
ylabel ( 'Radius  (ft)' 

zlabel ( 'Depth  [ft]') _ 
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APPENDIX  E.  SAMPLE  GEMGRID  CODES 


f  or  mat  =0  0  0  0  3 

<X  GRID  LI  NES> 

Dat  uml  ndex=l 
Dat  um=0. 

NCells=na  5 1  z  e  F I  r  s  t  Ce  I  I  =2  0 .  5 1  z  e  L  a  s  t  Ce  I  I  =2  0 .  Ratlo=na  Wl  dt  h=2  4  3  8  4.  #8  0  0ft 

<END  X  GRID  LI  NES> 

<Y  GRID  LI  NES> 

Dat  uml  ndex=l 
Dat  um=0. 

NCells=l  SlzeFlrstCell=na  SlzeLastCell=na  Ratlo=l.  Wldth=l.  #2D  Run 

#  NCells=na  5 1  z  e  F I  r  s  t  Ce  I  I  =2  0 .  5 1  z  e  L  a  s  t  Ce  I  I  =2  0 .  Ratlo=na  Wldth=9144.  #3  0  0ft 

<END  Y  GRID  LI  NES> 

<Z  GRID  LI  NES> 

Dat  uml  ndex=l 
Dat  um=-  1  6  0  0  2. 


NCel  1  s  =na 

SI  zeFI 

r  s  t  Ce  1  1 

=20. 

SI  zeLastCel  1 

=20. 

Rati  o=na 

Wl  dt  h=7  6  2. 

#25f  t 

clay 

NCe  1  1  s  =n  a 
water 

SI  zeFI 

r  s  t  Ce  1  1 

=20. 

SI  zeLastCel  1 

=20. 

Rati  o=na 

Wl  dt  h=1  5  2  4  0. 

#500f t 

NCel  1  s  =n  a 
<END  Z  GRI  D 

SI  zeFI 
LI  NES> 

r  s  t  Cel  1 

=20. 

SI  zeLastCel  1 

=20. 

Rati  o=na 

Wl  dt  h=7  6  2. 

#25f  t 

clay 

<OPTI  ONS> 

#P  I  ot  =Dy  s  ma  s  P 
<END  OPTI  ONS> 
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APPENDIX  F.  SAMPLE  PREGEMINI  CODES 


f  or  mat  =0  0  0  1  7 
<0PTI  0NS> 

Mode=start  #  PreGemini  mode  ( St  a  r  t  /  Rez  one/ Over  I  ay/ Repa  r  t  i  t  i  on ) 

StartTime=0.  #Set  initial  time  (REAL/precalc) 

Gravity=lg  #  Gravity,  positive  upward  (REAL/lg) 

<END  OPTI  ONS> 

<GRI  D> 

Coordi  nat  es =c y I i  n d r i  c a  I 
X Ce I  I  s  =1 2 1 9  d X  =0 .  x  Da t  u  m=0 

y  Ce I  I  s  =1  d y  =0 .  y  Da t  u  m=0 

z  Ce I  I  s  =8  3  8  d  z  =0 .  z  Da  t  u  m=0 

Gr  i  dFi  I  e  =  .  /  gr  i  d/  gr  i  d.  asc 

<END  GRI  D> 

<SUBGRI  DS> 

XSubGr  i  ds  =4  XPart  i  t i  on=aut  o 
YSubGr  i  ds  =1  YPart i  t i  on=aut  o 
ZSubGr  i  ds  =3  ZPart i  t i  on=aut  o 
<END  SUBGRI  DS> 

<BOUNDARY  CONDI  Tl  ONS> 

xmin=wall  xmax=free  #  wa  I  I  2  /  wa  I  I  /  f  r  e  e  n  g/ f  r  ee/ REAL 

#  y  mi  n  =wa  II  y  ma x  =f  r  ee 
z  mi  n  =wa  II  z  ma x  =f  r  ee 
<END  BOUNDARY  CONDI  Tl  ONS> 

<MATERI  ALS> 

Mat  er  I  al  I  D  =  he  solid 
Mat  er  i  al I  D  =  he’gas 
Ma  t  e  r  i  a  I  I  D=wa  t  e  r 
Mat  er I  al  I  D=ai  r 
Ma t  e  r I  a  I  I  D=c  I  a  y 
<END  MATERI  ALS> 

<BURN> 

Unburned=he  solid  Burned=he  gas  Time=0.  RefPt=(0.,  0.,  -  1  5  2  4.)  #50ft  z-deep  2D  Run 
<END  BURN> 

<HYDROSTATI  C  FI  ELD> 

pRef=1.0e-F6  zRef=0.  zMax=max  #  Ref  pressure,  ref  location,  zMax  (must  be 

1st  I  I  ne) 

Material=air  ei=eref  zMin=0.  #  top  state 
Mat  er  i  al  =wat  er  ei=eref  zMin=-  1  5  2  4  0.  #  2nd  state 

Ma  t  e  r  i  a  I  =c  I  a  y  ei=eref  zMin=-  1  6  0  0  2.  #  bottom  state 

<END  HYDROSTATI  C  FI  ELD> 


<1  Nl  Tl  AL  STATES>  #  State  matl  g  frac  rho  e  p  aO  u 

V  w 

StatelD  =  he  solid  Material=he  solid  Rho=rhoref  ei=eref  EOSvar=(0.)  # 

explosivestate 

St  at  el  D=wat  er  Mat  er  i  al  =wat  er  Rho=rhoref  ei=eref  EOSvar=(0.)  # 

water  state 

StatelD=air  Material=air  Rho=rhoref  ei=eref  EOSvar=(0.)  # 

at  mo  sphere  state 

StatelD=clay  Mat  er  i  al  =c  I  ay  Rho=rhoref  ei=eref  EOSvar=(0.) 

<END  I  Nl  Tl  AL  STATES> 

<FLOWFI  ELD> 

Opt  I  0 n  =h y  d  r  os  t  a t  i  c 

Option=ball  state=he  solid  mass=1  1  3  4  0.  RefPt=(0.,  0.,  -  1  5  2  4.)  #25lb  @50ft 

<END  FLOWFI  ELD> 

<TEXT  OUTPUT> 

imin=l  imax=l  jmin=l  jmax=l  kmin=l  kmax=l 

<END  TEXT  OUTPUT> 


File=hbx  Isolid.mtl 
Fi  I  e=hbx'l.  mtl 

"Fi  I  e=t  I  I  I  wat  er .  mt  I 
F  i  I  e  =a  i  r .  mt  I 
F  i  I  e  =c  I  a  y .  mt  I 


#  Coordi  nates:  CARTESIAN,  CYLINDRICAL,  or  SPHERICAL 

#  File  name  for  grid  file  specify  (none/STRING) 
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APPENDIX  G.  SAMPLE  GEMINI  CODES 


f  or  mat  =0  0  0  1  5 

<CASE> 

St  a  r  t  F  i  1  e 

=  ./pregemlnl/restart  000000  xxx.bln  #  Start  file  (use  "  xxx.bln 

ext  ens 1  on 

for  parallel 

runs) 

Tl  1 1  e=5  0  0  f  t 

# 

Title  (Limit  40  characters) 

Mode=f  1  ul 

d 

#  Run  mode  (Fluld/Coupled) 

<END  CASE> 

<TERMI  NATE> 

EndSt  ep=9  9  9  9  9  9 

#  Ter  ml  nat  1  on  step  (  1  NT/  none) 

EndTI  me=. 

15 

#  Term!  nation  time  (INT/none) 

dt  Ml  n  =1 .  e 

-12 

#  Terminate  If  step  size  Is  less  than  thi 

value 

(sec) 

<TRAP> 

#  x=6  0  8. 

44  I  =1  z 

=-  1  2  3.  8  6  Var=p 

Dl  f  =0.  01 

<END  TRAP> 

<END  TERMI  NATE> 

<1  NTEGRATI  0N> 

CFL=.  45 

#  CFL  safety  factor  (0.9  for  ID,  0.45  for 

2D  and  3D) 

CFLI  nl  t  =. 

05 

#  Initial  step  CFL  factor 

LI  ml  t  e r  =2 

#  LI  ml  ter:  0.0  -  2.0 

LagEqual  1  ze=sou 

#  Equalize  after  Lagrange  step 

Re  ma  p  Eq  u  a  1  1  z  e  =o  n 

#  Equalize  after  remap  step 

Pr  ot  ect  =1 

#  Protect  (  =0  off,  >0  on) 

<END  1  NTEGRATI  0N> 

<CELL  HI  ST0RY> 

1  =1 

j  =0.  k=- 

1  5  2  4. 

#charge  center  at  50ft  z-deep 

x=3  0  4  8. 

y=0. 

z=-  25. 

#at  X  =2 5cm  z-deep 

x=4  5  7  2. 

y=0. 

z=-  25. 

x=6  0  9  6. 

y=0. 

z=-  25. 

#pl  us  every  100ft  to  800ft 

x=7  6  2  0. 

y=0. 

z=-  25. 

x=9144. 

y=0. 

z=-  25. 

fexcept  at  150,  200,  250, 

3  0  0,  3  5  0, 

4  0  0ft 

x=1  0  6  6  8. 

y=0. 

z=-  25. 

x=12192. 

y=0. 

z=-  25. 

x=1  5  2  4  0. 

y=0. 

z=-  25. 

x=1  8  2  8  8. 

y=0. 

z=-  25. 

x=2  1  3  3  6. 

y=0. 

z=-  25. 

x=2  4  3  8  4. 

y=0. 

z=-  25. 

x=3  0  4  8. 

y=0. 

z=-  1  5  2  4. 

#a  t  X  =5  Of  t  z-deep 

x=4  5  7  2. 

y=0. 

z  =-  1  5  2  4. 

x=6  0  9  6. 

y=0. 

z=-  1  5  2  4. 

#pl  us  every  100ft  to  8  0  0ft 

x=7  6  2  0. 

y=0. 

z  =-  1  5  2  4. 

x=9144. 

y=0. 

z=-  1  5  2  4. 

x=1  0  6  6  8. 

y=0. 

z  =■  1  5  2  4. 

x=12192. 

y=0. 

z=-  1  5  2  4. 

x=1  5  2  4  0. 

y=0. 

z=-  1  5  2  4. 

x=1  8  2  8  8. 

y=0. 

z=-  1  5  2  4. 

x=2  1  3  3  6. 

y=0. 

z=-  1  5  2  4. 

x=2  4  3  8  4. 

y=0. 

z=-  1  5  2  4. 

x=3  0  4  8. 

y=0. 

z=-  2134. 

#a  t  X  =7  Of  t  z-deep 

x=4  5  7  2. 

y=0. 

z  =-  2134. 

x=6  0  9  6. 

y=0. 

z=-  2134. 

#pl  us  every  100ft  to  8  0  0ft 

x=7  6  2  0. 

y=0. 

z  =-  2134. 

x=9144. 

y=0. 

z=-  2134. 

x=1  0  6  6  8. 

y=0. 

z  =■  2134. 

x=12192. 

y=0. 

z=-  2134. 

x=1  5  2  4  0. 

y=0. 

z=-  2134. 

x=1  8  2  8  8. 

y=0. 

z=-  2134. 

x=2  1  3  3  6. 

y=0. 

z=-  2134. 

x=2  4  3  8  4. 

y=0. 

z=-  2134. 

<END  CELL  HI  ST0RY> 

<CONTOUR  PL0TS> 

St  epF  r  eq  = 

3  0  0  0 

#  Plot  step  frequency  (  1  NT/ n  o  n  e/ 1  a  b  1  e ) 

T 1  me  F  r  e  q  = 

0.  001 

#  Plot  time  frequency  (REAL/none/table) 

111 


<TARGET  PLOT  STEPS> 

3  0  0  0 

<END  TARGET  PLOT  STEPS> 

<TARGET  PLOT  Tl  MES> 

0.  0  0  0  1 
0.  001 
0.  025 
0.  030 
0.  035 
0.  040 
0.045 
0.  050 
0.  055 
0.  060 
0.  065 
0.  070 
0.  075 
0.  080 
0.  085 
0.  090 
0.  095 
0.  100 
0.  105 
0.  110 
0.  120 
0.  130 
0.140 
0.  150 

<END  TARGET  PLOT  Tl  MES> 

<END  CONTOUR  PLOTS> 


<RESTART> 

St  epF  r  eq  =1000 
SaveFreq=100 

yse  2  to  delete  every  other  file) 
<END  RESTART> 


#  Restart  file  output  step  frequency  (INT/none) 

#  Retention  Interval  (ex:  use  1  to  save  every  file, 


<TEXT  OUTPUT> 

St  epF  r  eq  =500 
Tl  meFreq=1000. 

Iniln=l  lniax=l  jniln=l  jniax=l  kmln^ 
<END  TEXT  OUTPUT> 


#  Output  step  frequency  (INT/none) 

#  Output  time  frequency  (INT/none) 
1  kmax=l 


#<COUPLED> 

#Un  I  t  s  =c  m-  g  -  s 
#Ba c  k P r  es  s  u  r  e  =0 .  e  +6 
#LoadFal  I  edSWI  =on 
#Bod  y  Or  I  g  I  n  =(  0 .  ,  0 .  , 
( c  m) 

#BodyXAxl  s=(  1.  ,  0., 
Fluid  coordinates  (cm) 
IBodyYAxI  s=(  0.  ,  1., 
Fluid  coordinates  (cm) 
#Bod y  ZAx  I  s  =(  0 .  ,  0 .  , 
Fluid  coordinates  (cm) 
#<END  COUPLED> 


#  Back  pressure  Inside  SWI  bodies  (d/cm''2) 

#  Load  failed  SWI  elements  (on/off) 

#  Location  of  Body  coordinate  system  In  Fluid  mesh 


s=(  1.  ,  0.  ,  0.  ) 
nates  (cm) 

# 

Vector 

along 

Body 

coordi 

nat  e 

system  x 

axis  i 

s=(  0.  ,  1.  ,  0.  ) 
nates  (cm) 

# 

Vector 

along 

Body 

coordi 

nat  e 

system  y 

axis  i 

s=(  0.  ,  0.  ,  1.  ) 
nates  (cm) 

# 

Vector 

along 

Body 

coordi 

nat  e 

system  z 

axis  i 

#<ELEMENT  HI  STORY> 

#1 

#2 

#3 

#<END  ELEMENT  HI  STORY> 


#<NODE  HI  STORY> 

#1 

#2 

#3 

#4 

#<END  NODE  HI  STORY> 


APPENDIX  H.  SAMPLE  GEMFIELD  CODES 


For  mat  =0  0  0  0  4 
<Ca  s  e  > 

DirOut=./  #  Subdirectory  for  output  files.  (Required.  Make  optional?) 

#  Use  "./"  for  current  working  directory. 

Series  =xyz  #  Series:  Output  file  identifier  (1-10  characters  or  ""). 

(  Opt i  0  n  a  I  ) 

F  i  I  e  F  0  r  ma  t  =Dy  s  ma  s  P  #  DysmasP  (latest  format),  DysmasP2010 

or  DysmasP  2  0  0  7 
# 

#  F  i  I  e  F  0  r  ma  t  =Tec  p  I  ot  Fi  I  eSpI  i  t  =Si  ngl  e  #  Tecpiot  (latest  format),  TecpIotlO, 

or  Tec  p I  ot  9 

#  F  i  I  e  F  0  r  ma  t  =Tec  p  I  ot  Fi  I  eSpI  i  t  =ByVar  #  Tecpiot  (latest  format),  TecpIotlO, 

or  Tec  p I  ot  9 

#  FileFormat=Tecplot  FileSplit=ByTime  #  Tecpiot  (latest  format),  TecpIotlO, 

or  Tec  pi  ot  9 

#  FileFormat=Tecplot  FileSplit=ByVarAndTime  #  Tecpiot  (latest  format),  TecpIotlO, 

or  Tec  p I  ot  9 

#  Fi  I  eFor  mat  =vt  kxml  F I  I  eTy  pe  =b  I  n  a  r  y  #  vtkxml  (ascii  or  binary)  or 

vtkiegacy  (ascii  only) 

#  Fi  I  eFor  mat  =vt  kxml  F  i  I  eTy  pe  =a  s  c  i  i  #  vtkxml  (ascii  or  binary)  or 

vtkiegacy  (ascii  only) 

#  F  i  I  e  F  0  r  ma  t  =vt  k  I  eg  a  c  y  #  vtkxml  (ascii  or  binary)  or 

vt  kl  egacy  (ascii  only) 

<E n d  Ca s e > 


#Note:  all  vars  are  optional 
<TI  ME> 

tBeg=-  1.0  0  0  tEnd=.25  #  tBeg,tEnd:  Time  window  for  output  (sec)  (before  scaling 

or  offset) 

#  Default:  Entire  s  i  mu  I  a  t  i  o  n  t  i  me 

tOffset=0.00  tScale=1.00  #  tOffset,  tScale:  Offset  and  scaling  for  time  (Default: 
t  Of  f  s  et  =0 ,  t  Sc  a  I  e  =1 ) 

#  t  <  =  =  t  Sc  a  I  e*  ( t  Of  f  s  et  -Ft ) 

MaxTi  meRec  or  ds  =5  0  0  0  0  #  Max  number  of  time  records  to  generate. 

#Note:  These  target  times  are  in  ADDITION  to  any  data  at  times  within  the  time 
window  (tBeg,  tEnd) 

#specified  above.  If  you  want  only  the  target  times  then  set  tBeg  and  tEnd  negative 
(or  very  large). 

<TI  ME  STEP  TARGETS> 

#1815 
#6  0  0  0 

<END  Tl  ME  STEP  TARGETS> 

<END  Tl  ME> 


<OPTI  ONS> 

S  h  0  wF  a  i  I  ed  =0  n  #  Show  failed  body  elements  (on  or  off,  default  is  on) 

#  (Use  "on"  to  preserve  correct  element  numbering  in  Tecpiot) 
Verbosity=l  #  Amount  of  screen  output  (Default=l) 

<END  OPTI  ONS> 


#  At  least  one  plot  variable  is  required 

#  Default  units  are  cgs 
<FLUI  D  VARI  ABLES> 

##These  vars  are  only  available  as  cell-averaged  quantities 

Var=p  Units=psi  #Units:  psi,  ksi.  Pa,  MPa,  bar  or  scaling  factor 

Var=u  Units=l.  #Units:  mis,  ft/s,  in/s  or  scaling  factor 

Var=v  Units=l.  #Units:  mis,  ft/s,  in/s  or  scaling  factor 

Var=w  Units=ft/s  #Units:  mis,  ft/s,  in/s  or  scaling  factor 

Var=ekin  #Units:  J,  MJ  or  scaling  factor 

##These  vars  are  available  as  c  el  I  -  aver  aged  (MatNum=0)  or  per  material  (MatNum>0) 

Var=r  MatNum=0  Units=1.0  #Units:  kg/m''3  or  scaling  factor 

Var=r  MatNum=l  Units=1.0  #Units:  kg/m""]  or  scaling  factor 

Var=r  MatNum=2  Units=1.0  #Units:  kg/m""]  or  scaling  factor 

Var=e  MatNum=0  Units=1.0  #Units:  J,  MJ  or  scaling  factor 

Var=e  MatNum=l  Units=1.0  #Units:  J,  MJ  or  scaling  factor 

Var=e  MatNum=2  Units=1.0  #Units:  J,  MJ  or  scaling  factor 

Va  r  =e i  n t  Mat  Nu  m=0  Un  i  t  s  =1 .  0 _ #Un  i  t  s :  J  ,  MJ  or  scaling  factor _ 
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Var  =ei  nt  Mat  Num=l 
Var  =ei  nt  Mat  Nuni=2 


Uni  ts=l.  0 
Un i  t  s  =1 .  0 


#Uni  ts: 
#Uni  ts: 


or  scaling  factor 
or  scaling  factor 


##These  vars  are  only  available  per  material  (MatNum>0) 
Var  =f  Mat  Nu m=l  #No  units 

Var  =f  Mat  Nu m=2  #No  units 

<END  FLUID  VARI  ABLES> 


<BODY  NODE  VARI  ABLES> 

#  Va  r  =u 

#  Va  r  =v 

#  Va  r  =w 

<END  BODY  NODE  VARI  ABLES> 


Uni  t  s  =1 . 
Uni  t  s  =1 . 


<BODY  ELEMENT  VARI  ABLES> 

#  Va  r  =p  Uni 

#  Var  =pPos  Uni 

#  Var  =pNeg  Uni 

<END  BODY  ELEMENT  VARI  ABLES> 


ml  i,  f  t  /  s , 
ml  s,  Ills, 
ml  s,  Ills, 


or  seal 
or  seal 
or  seal 


n  g  factor 
n  g  factor 
n  g  factor 


Pa,  MPa,  bar  or  seal 

Pa,  MPa,  bar  or  seal 

Pa,  MPa,  bar  or  seal 


ng  factor 
n  g  factor 
n  g  factor 


#Opt I  0  n  a  I  . 
<SUBDOMAI  N> 

I  Mi  n=l 
y  Mi  n  =0 . 
k  Mi  n  =1 
Offset=(0. 


If  left  unspecified,  the  entire  domain  will  be  output 


<END  SUBDOMAI  N> 


I  Max=9999999 
y  Ma  X  =0 . 
kMax=9999999 
0.  ,  0.  ) 


I  De  I  t  a  =1 
j  De  I  t  a  =1 
k  De  I  t  a  =1 


##################################################### 

#  The  following  is  a  Key  to  the  Plottable  Variables 
##################################################### 

#  Variable:  all  owed  Mat  Num  ( 0.  .  . 


# 

all  = 

all  point  variables 

NA 

# 

u  = 

velocity  in  r,x-dir 

0 

# 

V  = 

velocity  in  y,theta-dir 

0 

# 

w  = 

velocity  in  z-dir 

0 

# 

P  = 

pressure 

0 

# 

r  = 

dens  I  t  y 

a  n  y 

# 

e  = 

total  specific  energy 

a  n  y 

# 

el  nt  = 

internal  energy 

a  n  y 

# 

e  k  I  n  = 

kinetic  energy 

0 

# 

t  = 

t  e  mp  e  r  a  t  u  r  e 

0 

# 

f  = 

V  0 1  u  me  fraction 

>0 

# 

I  = 

I  mp  u  1  s  e  intensity  ( 1 1  me  in  sec) 

0 

# 

eos  #= 

EOS  variable  number  # 

0 

feature, 

# 

it  is  picked  automatically 

it- 

# 

St  at  = 

status  of  cell  ( ac  1 1  ve/  bl  oc  ked ) 

0 

# 

mat  = 

ma  t  e  r  I  a  1  id 

NA 

(for  his  only) 
(only  1  ma  t  e  r  I  a  I 


al  I  owed  to  have  this 


(currently  only  available  for  P-alpha  EOS) 
(for  field  only) 


########################################################## 

#  The  following  are  Scaling  Factors  for  some  common  units 
########################################################## 

#  Pressure 


1.  450377e- 005 

#ps  I 

1.  450377e- 008 

#ks  I 

0.  1 

#Pa 

1.  e-  007 

#MPa 

1.  e-  006 

#ba  r 

Force 

1.  e-  005 

#N 

1.  e-  008 

#kN 

2.  248089e- 006 

#1  bf 

#  Acceleration 


0.  01 

#m/  s  ''2 

0.  0328083 

#f  t/  0^2 

0.  3937008 

#1  n/  s  ''2 

Vel  oc  I  t  y 

0.  01 

#m/  s 

0.  0328083 

#f  t  /  s 

0.  3937008 

#1  n/  s 

Di  stance 

0.  01 

#m 

0.  0328083 

#f  t 

0.  3937008 

#1  n 

Vo  1  u  me 

1.  e-  006 

#m'' 

# 

3.5314676-005 

#f  t  "3 

# 

0.  06102374 

#i  n  "3 

# 

De  n  s i  t  y 

# 

1  0  0  0  #kq/m''3 

# 

Energy 

# 

1 .  e  -  0  0  7  #J 

# 

l.e-013  #MJ 

# 

Mass 

# 

0.001  #kg 

# 

1  mp  u  1  s  e  i  n  t  e  n  s  i  t  y 

# 

0.  1 

#Pa*s 

# 

1.  e-  007 

#MPa*s 

# 

1.  e-  006 

#ba  r  ♦  s 

# 

1.4  5  0  3  7  7  6  -  0  0  5 

#ps i  *6 

# 

1.4  5  0  3  7  7  6  -  0  0  2 

#ps i  ♦  ms 
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APPENDIX  I.  SAMPLE  GEMHIS  CODES 


For  mat  =0  0  0  0  4 
<Ca  s  e  > 

DirOut=./  #  Subdirectory  for  output  files.  (Required.  Make  optional?) 

#  Use  "./"  for  current  working  directory. 

Series=""  ISeries:  Output  file  identifier  (1-10  characters  or  ""). 

(  Opt i  0  n  a  I  ) 

F  i  I  e  F  0  r  ma  t  =d  y  s  ma  s  p 
or  DysmasP  2  0  0  7 
# 

#  F  i  I  e F 0 r  ma  t  =Tec  p I  ot 
for  mat),  TecpIotlO,  or 

#  F  i  I  e F 0 r  ma  t  =Tec  p I  ot 
for  mat),  TecpIotlO,  or 

#  F  i  I  e F 0 r  ma  t  =Tec  p I  ot 
or  Tec  pi  ot  9 

#  F  i  I  e F 0 r  ma  t  =Tec  p I  ot 
or  Tec  p I  ot  9 

#  F  i  I  e F 0 r  ma  t  =Tec  p I  ot 
or  Tec  pi  ot  9 

# 

#  FileFormat=plain 
<E n d  Ca s e > 


#  DysmasP  (latest  format),  DysmasP2010 


Fi  I  eSpI  i  t  Pt  =Si  ngl  e 
Tec  pi  ot  9 

Fi  I  eSpI  i  t  Pt  =Si  ngl  e 
Tec  pi  ot  9 

Fi  I  eSpI  i  t  Pt  =ByVar 

Fi  I  eSpI  i  t  Pt  =ByPt 

Fi  I  eSpI  i  t  Pt  =ByVar  AndPt 


FileSplitGbl=Single  #  Tecpiot  (latest 
Fi  I  eSpI  i  t  GbI  =ByVar  #  Tecpiot  (latest 

#  Tecpiot  (latest  format),  TecpIotlO, 

#  Tecpiot  (latest  format),  TecpIotlO, 

#  Tecpiot  (latest  format),  TecpIotlO, 


Fi  I  eSpI  I  t  Pt  =ByVar AndPt  # 


#Not  e:  all  vars  are  optional 
<Ti  me> 

tBeg=-  1.0  0  0  tEnd=.15 
or  offset) 

t  Of  f  s  et  =0.  00  t  Sc  a  I  e  =1 .  00 
t  Of  f  s  et  =0 ,  t  Sc  a  I  e  =1 ) 

Ma xTi  me Rec  0 r  d s  =5  0  0  0  0  0 
(default=5  0  0  0  0  0) 

RemoveOver  I  a  p=on 

overlap  (on  or  off,  default  is 

<End  Ti  me> 


#  tBeg,tEnd:  Time  window  for  output  (sec)  (before  scaling 

#  Default:  Entire  s  I  mu  I  a  1 1  o  n  1 1  me 

#  tOffset,  tScale:  Offset  and  scaling  for  time  (Default: 

#  t  <  =  =  t  Sc  a  I  e*  ( t  Of  f  s  et  -Ft ) 

#  Set  max  num  of  time  records  (used  for  array  allocation) 

#  Flag  for  removing  oldest  data  if  there  is  a  time 
on) 


#  Default  units  are 
<Point  Variables> 
##Th es e  vars  are  o n I 
Va  r  =p 
Va  r  =u 
Va  r  =v 
Va  r  =w 
Va  r  =e k I  n 
Va  r  =1 


cgs 


y  available  as 
Uni  ts=l. 

Uni  ts=l. 

Uni  ts=l. 

Uni  ts=l. 


##These 

vars 

are  avail 

able 

as  cell 

Va  r  =r 

Mat  Nu  m=0 

Un 

ts=l. 

0 

Va  r  =r 

Mat  Nu  m=l 

Un 

ts=l. 

0 

Va  r  =r 

Mat  Nu  m=2 

Un 

ts=l. 

0 

Va  r  =e 

Mat  Nu  m=0 

Un 

ts=l. 

0 

Va  r  =e 

Mat  Nu  m=l 

Un 

ts=l. 

0 

Va  r  =e 

Mat  Nu  m=2 

Un 

ts=l. 

0 

#Va  r  =e 

Mat  Nu  m=3 

Uni  t  s  =1 

.  0 

Va  r  =e i 

nt 

Mat  Nu  m=0 

Un 

ts=l. 

0 

Va  r  =e I 

nt 

Mat  Nu  m=l 

Un 

ts=l. 

0 

Var  =ei 

nt 

Mat  Nu  m=2 

Un 

ts=l. 

0 

cell-averaged  quantities 

#Units:  psi,  ksi.  Pa,  MPa,  bar  or  scaling  factor 
#Units:  mis,  ft/s,  in/s  or  scaling  factor 

#Units:  mis,  ft/s,  in/s  or  scaling  factor 

#Units:  mis,  ft/s,  in/s  or  scaling  factor 

#Un  i  t  s :  J  ,  MJ  or  scaling  factor 
#No  units 

averaged  (MatNum=0)  or  per  material  (MatNum>0) 


#Uni  t  s 

kg/  m''3 

or 

seal 

I  ng 

factor 

#Uni  t  s 

kg/  m''3 

or 

seal 

I  nc 

factor 

#Uni  t  s 

kg/  m''3 

or 

seal 

I  nc 

factor 

#Uni  t  s 

J  .  MJ 

or 

seal  i 

ng 

factor 

#Uni  t  s 

J  .  M] 

or 

seal  i 

ng 

factor 

#Uni  t  s 

1  ,  M] 

or 

seal  i 

ng 

factor 

#Uni  ts:  1  ,  Ml 

or 

seal 

I  n  g 

factor 

#Uni  t  s 

J  .  M] 

or 

seal  i 

ng 

factor 

#Uni  t  s 

J  .  M] 

or 

seal  I 

ng 

factor 

#Uni  t  s 

J  .  M] 

or 

seal  I 

ng 

factor 

##These  vars  are  only  available 
Var=f  MatNum=l 

Var=f  MatNum=2 

<End  Point  Variables> 


per  ma  t  e  r  i  a  I  (Mat  Nu  m>0 ) 
#No  units 
#No  units 


#################################### 
#The  remaining  sections  are  optional 
#################################### 
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<Opt I  ons  > 

Ch  a  r  g  eSc  a  I  e  =  1 .  #  Scaling  for  specie  (Energyint,  Energy,  Mass,  Volume,  Radius) 

( def  aul  t  =1.  ) 

Verbosity=l  #  IVerbose:  Screen  output  (0=no,  l=some,  2=lots) 

pRef  I  mp=- 5.  OOe?  #  Impulse  intensity  ref  pressure  ( d  y  n  es  /  c  m''2 )  (if  unset  then  use 

value  from  1st  record  in  file) 

<E n d  Opt i  0 n s  > 

<Global  Variables> 

#Note:  MatNum=0  is  for  total  of  all  materials 

Var=radius  MatNum=l  Units=l.  #Units:  m,  ft,  in  or  scaling  factor 

Var=volcav  MatNum=3  Units=l.  #Cavitation  volume  in  water 

#for  body/block  fixed  mesh  option 

#  Var=BodyCG-X 

#  Var=BodyCG-Y 

#  Var=BodyCG-Z 

#  Var =BodyVel  -  X 

#  Var =BodyVel - Y 

#  Var =BodyVel  ■  Z 

#  Var=BodyF-X 

#  Var=BodyF-Y 

#  Var=BodyF-Z 

<End  Global  Variables> 


<Body  Global  Var  I  abl  es  > 

Va  r  =f  0 r  c  e-  X 
Va  r  =f  0 r  c  e-  y 
Va  r  =f  or  c  e-  z 
Va  r  =f  0  r  c  e 
#Va  r  =1  mp u  I  s  e 

<End  Body  Global  Variables> 

<Body  Element  Variables> 

Va  r  =p 
#  Var=i 

<End  Body  Element  Variables> 

<Body  Node  Variables> 

Va  r  =f  0 r  c  e-  X 
Va  r  =f  0 r  c  e-  y 
Va  r  =f  0 r  c  e-  z 
Va  r  =f  0  r  c  e 
Var  =u 
Va  r  =v 
Va  r  =w 
Var  =vel 
Var  =x 
Va  r  =y 
Va  r  =z 

<End  Body  Node  Variables> 


##################################################### 

#  The  following  is  a  Key  to  the  Global  Variables 
##################################################### 

#  Variable:  all  owed  Mat  Num  ( 0.  . 

#all=  all  global  variables  NA 

#dt=Geminitimestepsize  NA 

#  ma  s  s  =  ma  s  s  any 

#energy=total  energy  any 

#radius=equivalent  radius  >0 

#  vol  =  vol  u  me  >0 

#volcav=cavitatedvolume  >0 

#BodyVel-[xyz]  NA  (forbo 

#BodyCG-[xyz]  NA  (forbo 

#  BodyF-  [xyz]  NA  (forbo 


(for  body/block  f 
(for  body/block  f 
(for  body/block  f 


xed  mesh  opt 
xed  mesh  opt 
xed  mesh  opt 


##################################################### 
#  The  following  is  a  Key  to  the  Point  Variables 
##################################################### 


#  Va  r I  a  b 


a  1 

1 

poi  nt 

var 

I  abl  es 

vel 

oc 

ity  I 

n  r , 

X  -  d  I  r 

vel 

oc 

ity  I 

n  y- 

t  het  a- 

vel 

oc 

ity  I 

n  z  - 

di  r 

al  I  owed  Mat  Num  ( 0.  .  .  NMat ) 
NA 
0 
0 
0 


# 

p  = 

pressure 

0 

# 

r  = 

dens  i  t  y 

a  n  y 

# 

e  = 

t  ot  a i  specific  energy 

a  n  y 

# 

ei  nt  = 

internai  energy 

a  n  y 

# 

e  k  i  n  = 

kinetic  energy 

0 

# 

t  = 

t  e  mp  e  r  a  t  u  r  e 

0 

# 

f  = 

V  0  i  u  me  fraction 

>0 

#  eos  # 
feature 

# 

#  St  at  = 

#  mat  = 


i  mp  u  I  s  e  i  I 
=  EOS  va  r i i 
it  is  p  i  I 

status  of 
ma  t  e  r  i  a  i  i 


n  t  e  n  s  i  t  y  ( t  i  me  in  sec) 
abi  e  number  # 
eked  automaticaiiy 

cei  i  ( act  i  ve/  bi  ocked) 


(for  his  oniy) 

(oniy  1  materiai  aiiowed  to  have  this 

(currentiy  oniy  avaiiabie  for  P-aipha  EOS) 
(for  fieid  oniy) 


##################################################### 

#  The  foiiowing  is  a  Key  to  the  Body  Giobai  Variabies 

##################################################### 
#force-x=Totai  body  force  in  r,x-dir  (sum  of 

#force-y=Tota(  body  force  in  y,theta-dir  (sum  of 

#force-z=Tota(  body  force  in  z-dir  (sumof 

#  force  =Totai  body  force  magnitude  (sum  of 

#  i  mp  u  i  s  e  =  i  mp  u  i  s  e  ( t  i  me  in  sec) 

##################################################### 

#  The  foiiowing  is  a  Key  to  the  Body  Eiement  Variabies 
##################################################### 

#  p  =  Eiement  pressure  (net) 

#  i  =  i  mp  u  i  s  e  intensity  ( t  i  me  in  sec) 

#  ppos  =  Eiement  pressure  on  "positive"  side  (normai 
interface  eiement 

#  pneg  =  Eiement  pressure  on  "negative"  side  (normai 
coupiing  interface  eiement 

##################################################### 

#  The  foiiowing  is  a  Key  to  the  Body  Node  Variabies 
##################################################### 


node  forces) 
node  forces) 
node  forces) 
node  forces) 


points  toward  "  ey 
points  away  from 


of  coupi  i  ng 


#  f  0 r  c  e-  X 

#  f  or  c  e-  y 

#  f  0 r  c  e-  z 

#  force 

#  u  =  N 

#  V  =  N 

#  w  =  N 

#  V  e  i  =  N 

#  X  =  N 

#  y  =  N 

#  z  =  N 


Node  force 
Node  force 
Node  force 


n  r ,  X  -  d  i  r 
n  y,  t  het  a- di  r 
n  z-dir 


Node  force  magnitude 


Node  vei  oc 
Node  vei  oc 
Node  vei  oc 
Node  vei  oc 
Node  posit 
Node  posit 


Node  pos 


t  y  in  r , X  -  d  i  r 
ty  in  y,theta-dir 
t  y  in  z-dir 
t  y  ma  g  n  i  t  u  d  e 
on  ( r , X - d i  r ) 
on  ( y,  t  het  a- di  r ) 
on  (z-dir) 


########################################################## 

#  The  foiiowing  are  Scaiing  Factors  for  some  common  units 
########################################################## 

#  Pressure 

#  1.  4  5  0  3  7  7  e-  0  0  5  #psi 

#  1.  4  5  0  3  7  7  e-  0  0  8  #ksi 

#  0.1  #Pa 

#  l.e-  0  0  7  #MPa 

#  l.e-  0  0  6  #bar 

#  Force 

#  l.e-  0  0  5  #N 

#  l.e-  0  0  8  #kN 

#  2.  2  4  8  0  8  9  e-  0  0  6  #i  bf 

#  Acceieration 

#  0.01  #m/s''2 

#  0.0  3  2  8  0  8  3  #ft/s''2 

#  0.3  9  3  7  0  0  8  #in/s''2 

#  Vei  oc i  t  y 

#  0.01  #m/s 

#  0.0  3  2  8  0  8  3  #ft/s 

#  0.3  9  3  7  0  0  8  #in/s 

#  Distance 

#  0.01  #m 

#  0.0  3  2  8  0  8  3  #ft 

#  0.3  9  3  7  0  0  8  #i  n 

#  Vo i  u  me 

#  l.e-  0  0  6  #m''3 

#  3.  5  3  1  4  6  7  e-  0  0  5  #ft''3 

#  0.0  6  1  0  2  3  7  4  #in''3 

#  Density 

#  1  0  0  0  Ikg/m"! 


APPENDIX  J.  SAMPLE  MATERIAL  FILE 


f  or  mat  =0  0  0  0  2 
<0PTI  0NS> 

title="sand  eos  using  miegrunisen  eos" 
eosType=mi  eg 
<END  OPTI  0NS> 

<REFERENCE> 
r  hoRef  =2.  023 
e i  Ref  =0 . 
cRef  =2  6  0  0  0  0.  0 
<END  REFERENCE> 

<LI  Ml  TS> 

r  hoMi  n=l.  e-  02 
ei  Mi  n=-  9.  9e+99 
p Mi  n  =2 .  e  +4 
<END  LI  Ml  TS> 

<E0S  VARS> 
g  a  mma  0  =  .  9  7 
5=1.  86 
rho0=2.  023 
ei  0=0.  0 

p0=0.  0 

c  0  2  =4.  0  4  8  e+10 
<END  EOS  VARS> 
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